使用 matlab 进行傅里叶分析和滤波-技术圈

作者：盛夏的猫宁

来源：SegmentFault 思否社区

傅里叶分析

公式法

下例是将振幅为1的5Hz正弦波和振幅为0.5的10Hz正弦波相加之后进行傅里叶分析。

clear allN=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;endfor 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);    endend
%绘图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)。x是原始信号，y是变换之后的信号。y与x有相同的长度。
y=fft(x,N)。x，y定义如上，N是正整数，表示进行N点快速傅里叶变换。如果x长度小于N，则对x补零，使之与N相等；反则，则对x进行截取。

对应的逆变换有两种，分别为x=ifft(y)和x=ifft(y.N)。

一般而言，N点fft的结果y，在处对应的频率为最高采样率的一半，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);    endend
%逆变换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，则）。可以直接给定，亦可以根据参数由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，则。可以直接给定，亦可以根据参数由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)
其他滤波器，如椭圆滤波器等，使用方式类似，只是函数名称不同。

-END-
点击左下角阅读原文，到 SegmentFault 思否社区 和文章作者展开更多互动和交流。