【离散数据处理】快速傅里叶变换与滤波

快速傅里叶变换（FFT）与滤波是数字信号处理中很成熟的方法，相关教材很多，记录一下常用方法。

FFT

例如有测量数据集 $\{(t_i,y_i)|i=1,\dots ,n\}$，由以下函数生成：
$$\begin{array}{rl}{y}_1& =2\cdot sin(2\pi f_1\cdot t)\\ y_2& =sin(2\pi f_2\cdot t+\phi )\\ y& =y_1+y_2\end{array}$$
$f_1=1Hz、f_2=5Hz$， 采样频率 $f_s=100Hz$，FFT结果如下：

5个周期数据	幅频曲线
相频曲线(phi=-90°)	相频曲线(phi=90°)

• 幅频曲线看到有$f_1=1Hz、f_2=5Hz$的两个信号。
• 某些领域较少关注FFT的相频曲线，我也没完全搞懂。由于函数angle返回的是余弦的相位，而 $cos0°=sin(90°)$，达到相同的值，余弦比正弦提前1/4个周期，也就是 $sint=cos(t-90°)$，理论上返回的是 $\phi -90°$，然而当 $\phi =90°$ 时并未返回0°，感兴趣的可以试试DFT。

滤波

通过巴特沃斯滤波器来实现低通滤波器。

滤波器幅频特性

滤波前后比较

• 经过滤波器后，高频率信号被去除，但是滤波后的信号往后延迟了1/4个周期左右，相移量受滤波器参数控制，这是巴特沃斯滤波器的不足，有空探讨一下其它的方法。

详细代码

%
% Auther:  sanfan66
%
clc
clear all
close all

fs=100;
t=0:1/fs:5;
f1=1;
f2=5;
y1=2*sin(2*pi*f1*t);
y2=sin(2*pi*f2*t-pi/2);
y=y1+y2;

figure
plot(t,y1,t,y2,t,y)
legend('y1','y2','y')
xlabel('t')
ylabel('y')
xlim([t(1) t(end)]);
grid on;grid minor;
saveas(gcf, '1 y', 'png')

%% fft
num=length(t);
f=((0:num-1)/num-1/2)*fs;%双边谱频率
% f=linspace(-fs/2,fs/2,num);方法2
fftY=fft(y);
fftYShift=fftshift(fftY);
mag=abs(fftYShift)/(num/2);
phase=angle(fftYShift);%-pi*f/num

figure
plot(f,mag,'-o');
xlabel('f[hz]')
ylabel('magnitude')
xlim([0 fs/2])
grid on;grid minor;
saveas(gcf, '2 f-mag', 'png')
figure
plot(f,phase*180/pi,'-o');
xlabel('f[hz]')
ylabel('phase')
xlim([0 fs/2])
grid on;grid minor;
saveas(gcf, '3 f-phase', 'png')

%% 滤波
wp=2;%通带截止频率
ws=4;%阻带截止频率
rp=3;%通带最大衰减
rs=30;%阻带最小衰减

[n,wn]=buttord(wp/(fs/2), ws/(fs/2), rp, rs);
[b,a]=butter(n,wn,'low');
[h,w]=freqz(b,a,1024,fs);
filY=filter(b,a,y);

figure
plot(w,abs(h));
xlabel('f[hz]')
ylabel('滤波器特性')
grid on;grid minor;
saveas(gcf, '5', 'png')
figure
plot(t,y,t,y1,t,filY);
xlabel('t')
ylabel('y')
legend('y','y1','fil')
grid on;grid minor;
saveas(gcf, '6', 'png')