matlab運用fft進行簡單的濾波

******該程序共分為四個部分，繪制原始時域信號圖；進行傅里葉變換並繪制頻譜圖；濾波並繪制頻譜圖；進行傅里葉反變換，繪制濾波后時域信號圖******

clear all;
dt=0.02;N=512;
n=1:N-1;t=n*dt;                                         %時間序列
f=n/(N*dt);                                             %頻率序列
f1=3;f2=10;                                             %已知原始信號的頻率成分
x=0.5*sin(2*pi*3*n*dt)+cos(2*pi*10*n*dt);               %原始信號
subplot(2,2,1);plot(t,x);                               %時域信號圖
title('時域信號');xlabel('時間/s');


y=fft(x);                                               %對時域信號進行FFT變換
xlim([0 12]);ylim([-1.5 1.5]);
subplot(2,2,2);plot(f,abs(y)*2/N);                      %繪制原始信號的振幅圖
xlabel('頻率/Hz');ylabel('振幅');
xlim([0 50]);title('原始振幅譜');
ylim([0 0.8]);


f1=8;f2=15;                                             %要濾除頻率的上下限
yy=zeros(1,length(y));                                  %設置與y相同的元素數組
for m=0:N-1                                             %將頻率落在該頻率范圍內及其大於奈奎斯特的波濾除
    %小於奈奎斯特頻率的濾波范圍
    if (m/(N*dt)>f1&m/(N*dt)<f2)|(m/(N*dt)>(1/dt-f2)&m/(N*dt)<(1/dt-f1));
        %大於奈奎斯特頻率的濾波范圍
        %1/dt為一個頻率周期
        yy(m+1)=0;                                      %將落在此頻率范圍內的振幅設置為0
    else
        if m<511;                                       %確定數組y(m+1)不溢出
        yy(m+1)=y(m+1);                                 %其余頻率范圍內的信號保持不變
        end
    end
end
subplot(2,2,4);plot(f,abs(yy)*2/N);                     %繪制濾波后的振幅譜
xlim([0 50]);ylim([0 0.5]);
xlabel('頻率/Hz');ylabel('振幅');
gstext=sprintf('自%4.1f-%4.1fHz的頻率被濾除',f1,f2);
%將濾波范圍顯示作為標題
title(gstext);


subplot(2,2,3);plot(t,real(ifft(yy)));
%繪制濾波后的數據運用ifft變換回時域並繪圖
title('通過IFFT回到時域');
xlabel('時間/s');
ylim([-0.6 0.6]);xlim([0 12]);