Fs = 1000; % 采样频率
N = 1024; % 采样点数 ,即实际中一次FFT变换所使用的点数
n = 0:N-1; % 采样序列,为plot绘图用的序列,其对应点表示的真实频率为 f = f = n*(Fs/N),Fs/N,是频率分辨率,其中0 <= N < max(N)/2-1
f = (-N/2:N/2-1)*Fs/N;%真实频率序列,真实频率从(-N/2~(N/2)-1)*(Fs/N)
t = 0:1/Fs:(N/Fs)-1/Fs; % 时间序列
x = 2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*200*t+pi*90/180) ; %原始信号序列
subplot(2,2,1);%绘图 2行2列的第个1框
plot(t,x);%plot(x,y) 以x 元素为横坐标值,y 元素为纵坐标值绘制曲线。
grid on;
xlabel('t');
ylabel('y');
title('原始信号');
y = fft(x,N); %对原始信号做N点FFT变换
z = fftshift(y);%对fft结果进行转换,将频率移动的正确的位置
M = abs(z);
subplot(2,2,2);%绘图 2行1列的第2个框
plot(f,M);
grid on;
xlabel('Hz');
ylabel('幅度');%FFT变换后幅度值为原始信号对应频率的幅值*Fs
title('原始信号FFT');
b = fir1(28,0.25);%更多用法查fir1函数说明。fir1(n,Wn),默认选用哈明窗,n是滤波器阶数,Wn是一个归一化频率,在滤波器设计中是用fs/2进行归一。在fir1滤波器设计时采用的是归一化频率。实际采样频率为fs,实际的截止频率为fc,设归一化截止频率为fcm,fcm=fc/(fs/2)。
fir_output = filter(b,1,x);
subplot(2,2,3);%绘图 2行1列的第3个框
plot(t,fir_output);
grid on;
xlabel('t');
ylabel('y');
title('FIR滤波后信号(截止频率为 125 Hz)');
y2 = fft(fir_output,N); %对原始信号做N点FFT变换
z2 = fftshift(y2);%对fft结果进行转换,将频率移动的正确的位置
M2 = abs(z2);
subplot(2,2,4);%绘图 2行1列的第4个框
plot(f,M2);
grid on;
xlabel('Hz');
ylabel('幅度');
title('FIR滤波后信号FFT');