内容简介:下例 是将振幅为1的5Hz正弦波和振幅为0.5的10Hz正弦波相加之后进行傅里叶分析。运行结果如下所示:
傅里叶分析
公式法
下例 是将振幅为1的5Hz正弦波和振幅为0.5的10Hz正弦波相加之后进行傅里叶分析。
clear all N=512; dt=0.02; n=0:N-1; t=n*dt; x=sin(2*pi*5*t)+0.5*sin(2*pi*10*t);%生成和信号 %傅里叶变换 m = floor(N/2)+1; a=zeros(1,m); b=zeros(1,m); for k=0:m-1 for ii=0:N-1 a(k+1) = a(k+1)+2/N*x(ii+1)*cos(2*pi*k*ii/N); b(k+1) = b(k+1)+2/N*x(ii+1)*sin(2*pi*k*ii/N); end c(k+1)=sqrt(a(k+1).^2+b(k+1).^2); end %傅里叶逆变换 if(mod(N,2) ~=1) a(m)=a(m)/2; end for ii=0:N-1 xx(ii+1)=a(1)/2; for k=1:m-1 xx(ii+1)=xx(ii+1)+a(k+1)*cos(2*pi*k*ii/N)+b(k+1)*sin(2*pi*k*ii/N); end end %绘图 subplot(3,1,1),plot(t,x,'LineWidth',2);title('原始信号'),xlabel('时间/s'); subplot(3,1,2),plot((0:m-1)/(N*dt),c,'LineWidth',2);title('傅里叶变换'),xlabel('频率/Hz'); subplot(3,1,3),plot((0:N-1)*dt,xx,'LineWidth',2);title('合成信号'),xlabel('时间/s');
运行结果如下所示:
快速傅里叶
matlab中的快速傅里叶有两种调用形式:
y=fft(x) y=fft(x,N)
对应的逆变换有两种,分别为 x=ifft(y)
和 x=ifft(y.N)
。
一般而言,N点fft的结果y,在$n=N/2+1$处对应的频率为最高采样率的一半,y的后一半与前一半对称。
下例 是将振幅为1的5Hz正弦波和振幅为0.5的10Hz正弦波相加之后进行傅里叶分析。
clc;clear; fs=30; N=512; n=0:N-1; t=n/fs; x=sin(2*pi*5*t)+0.5*sin(2*pi*10*t); %快速傅里叶 y=fft(x,N); mag = abs(y);%振幅 ang = angle(y);%相位 f = n*fs/N;%横轴各点对应的频率值 %逆变换 xx = ifft(y); xx = real(xx);%计算误差使得xx可能是复数,对其取实部得到信号 tt = [0:length(xx)-1]/fs;%横轴各点对应的时间
结果图省略。
滤波
利用快速傅里叶简单滤波
下例是将振幅为1的5Hz正弦波和振幅为0.5的10Hz正弦波相加之后,滤除8Hz以上的信号。
clc;clear; fs=30;%采样率 N=256; n=0:N-1; t=n/fs; dt=1/fs; x=sin(2*pi*5*t)+0.5*sin(2*pi*10*t); %快速傅里叶 y=fft(x,N); mag = abs(y)*2/N;%振幅 ang = angle(y);%相位 f = n*fs/N;%横轴各点对应的频率值 %滤波 nn=length(y); yy = zeros(1,nn); for m=0:nn-1 if(m/(nn*dt)>8)&(m/(N*dt)<(1/dt-8)) yy(m+1)=0; else yy(m+1)=y(m+1); end end %逆变换 xx = ifft(yy); xx = real(xx);%计算误差使得xx可能是复数,对其取实部得到信号 tt = [0:length(xx)-1]/fs;%横轴各点对应的时间 %绘图 subplot(2,1,1),plot(t,x,'LineWidth',2);title('原始信号'),xlabel('时间/s'); subplot(2,1,2),plot(tt,xx,'LineWidth',2);title('滤波后的信号'),xlabel('时间/s');
结果如下图
几个简单的函数
-
xi=filtic(B,A,ys)
。B、A分别为系统 z变换 后的传递函数的分子和分母多项式的系数向量,ys是系统的初始输出状态,xi为对应的初始条件下输入序列。 -
yn0=filter(B,A,xn)
。B、A定义同上,xn是系统的输入信号,yn0为系统的零状态响应。 -
yn=filter(B,A,xn,xi)
。B、A、xn、xi定义同上,yn为系统全响应。
模拟滤波器
以巴特沃斯低通滤波器为例,说明调用方法。
[btt1,ctt1] = butter(N,wn,'s');%1.调用函数生成滤波器系数 H = [tf(btt1,ctt1)];%滤波器的传递函数 t = (0:n-1)./fs;%时域信号横轴的坐标,n为长度,fs为采样率 s1 = lsim(H,a1,t);%2.滤波
说明:
-
[btt1,ctt1] = butter(N,wn,'s');
。N是滤波器的阶数,wn是截止频率(是弧度值,如果截止频率要求为500Hz,则$wn=500*2*\pi$)。可以直接给定,亦可以根据参数由buttord`函数计算得到。's'表示模拟滤波器。btt1、ctt1分别表示滤波器在 拉普拉斯域 中传递函数的分子、分母多项式的系数。 -
s1 = lsim(H,a1,t)
。H是模拟滤波器的传递函数,a1表示待滤波信号,t是信号的横坐标,s1是滤波后的信号。
其他说明:
- 这里仅以低通滤波器为例,其他巴特沃斯滤波器如高通、带通、带阻调用方式类似,只是函数
butter
的参数略有不同,请参看matlab关于butter函数的介绍。(在matlab中执行help butter
) - 其他滤波器,如椭圆滤波器等,使用方式类似,只是函数名称不同。
数字滤波器
以巴特沃斯低通滤波器为例,说明调用方法。
%方式一:直接设计 [btt,ctt] = butter(N,wn);%1.生成数字滤波器 Signal_Filter=filter(btt,ctt,a1);%2.滤波 %方式二:模拟滤波器转数字滤波器 [btt1,ctt1] = butter(N,Wn,'s'); [btt1,ctt1]=bilinear(btt1,ctt1,Fs);%3.模拟滤波器转数字滤波器 Signal_Filter=filter(btt,ctt,a1);
说明:
-
[btt,ctt] = butter(N,wn)
。N是滤波器阶数,wn是相对截止频率,比如最高采样率为Fs,要求的截止频率为fs,则$wn=fs/Fs$ 。可以直接给定,亦可以根据参数由buttord
函数计算得到。注意,这里 没有 参数's'。btt、ctt分别表示滤波器在 z域 中传递函数的分子、分母多项式的系数。 -
Signal_Filter=filter(btt,ctt,a1)
。btt、ctt与之前定义相同,a1是待滤波信号,Signal_Filter是滤波之后的信号。 -
[btt1,ctt1]=bilinear(btt1,ctt1,Fs)
。是使用双线性法将模拟滤波器在 拉普拉斯域 中的系数转换成数字滤波器在 z域 中的系数,Fs是采样率。
其他说明:
- 这里仅以低通滤波器为例,其他巴特沃斯滤波器如高通、带通、带阻调用方式类似,只是函数
butter
的参数略有不同,请参看matlab关于butter函数的介绍。(在matlab中执行help butter
) - 其他滤波器,如椭圆滤波器等,使用方式类似,只是函数名称不同。
以上就是本文的全部内容,希望对大家的学习有所帮助,也希望大家多多支持 码农网
猜你喜欢:- openCV之中值滤波&均值滤波(及代码实现)
- 粒子滤波Matlab示例
- 粒子滤波Matlab示例
- OpenCV 线性滤波
- 【信号与系统】05 - 滤波、采样和通信
- C++ 实现 CVPR 2019 side window 中值滤波
本站部分资源来源于网络,本站转载出于传递更多信息之目的,版权归原作者或者来源机构所有,如转载稿涉及版权问题,请联系我们。