前言
对数字信号处理课程中的一些理论产生兴趣，通过Matlab进行一些实操，实现信号滤波的基本原理为傅立叶变换，对信号进行傅立叶变换后，则可得到频谱分布，依次便可进行滤波，包括高通、低通、带通等等。
基础函数
FFT（快速傅立叶变换）
IFFT（逆傅立叶变换）
X=FFT(x)；
X=FFT(x,N);
x=IFFT(X);
x=IFFT(X,N);
功能实现
创建数个间歇波叠加形成的信号，通过滤波处理，并显示处理前后的图像。
低通滤波部分的实现
clc
clear
close all;
%input signal-------------------------
Fs =1000;
T=1/Fs;
L=1000;
t=(0:L-1)*T;
y=0.7*sin(2*pi*10*t)+2*sin(2*pi*120*t);
figure(1)
subplot(3,2,1);
plot(t,y)
grid on
%side ft----------------------------------
NFFT=L;
Y=fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
subplot(3,2,2);
plot(f,2*abs(Y(1:NFFT/2+1)));
grid on;
%low filter---------------------------
fn =1002;
fp=60;
fs =100;
Rp=2;
Rs=30
Wp=fp/(fn/2);
Ws=fs/(fn/2);
[n,Wn]=buttord(Wp,Ws,Rp,Rs);
[b,a]=butter(n,Wn);
[H,F]=freqz(b,a,NFFT/2+1,1002);
% subplot(3,3,3);
% plot(F,20*log10(abs(H)));
grid on;
% subplot(3,3,4);
pha=angle(H)*180/pi;
% plot(F,pha);
grid on;
A = 2*abs(Y(1:NFFT/2+1)).*(abs(H)');
subplot(3,2,3);
plot(f,A);
grid on;
% %ifft-------------------------------
for(i =1:L-(NFFT/2+1))
A(1,((NFFT/2+1)+i))=A(1,((NFFT/2+1)-i));
end
z=A.*exp(1i*angle(Y));
z2=real(ifft(z));
subplot(3,2,4);
y=0.7*sin(2*pi*10*t);
plot(t,y);
axis([0,1,-2,2]);
% subplot(3,2,7);
% plot(t,z2*NFFT/2);
%high filter--------------------------
fn=1002;
fp=60;
fs=30;
Rp=3;
Rs=20;
Wp=fp/(fn/2);
Ws=fs/(fn/2);
[n,Wn]=buttord(Wp,Ws,Rp,Rs);
[b,a]=butter(n,Wn,'high');
[H,F]=freqz(b,a,NFFT/2+1,1002);
pha=angle(H)*180/pi;
A=(2*abs(Y(1:NFFT/2+1))).*(abs(H)');
subplot(3,2,5);
plot(f,A);
grid on;
%ifft-----------------------------
for i=1:L-(NFFT/2+1)
A(1,((NFFT/2+1)+i))=A(1,((NFFT/2+1)-i));
end
z=A.*exp(1i*angle(Y));
z2=real(ifft(z));
subplot(3,2,6);
y = sin(2*pi*120*t);
plot(t,y);
axis([0,1,-2,2]);
带通部分的实现
clc
clear
close all
%????????
Fs = 1000;
T=1/Fs;
L=1000;
t=(0:L-1)*T;
y=0.7*cos(2*pi*100*t)+sin(2*pi*200*t)+0.4*sin(2*pi*40*t)+5*cos(2*pi*400*t);
subplot(3,2,1)
plot(t,y)
grid on
NFFT=L;
Y=fft(y,NFFT)/L;
f=Fs/2*linspace(0,1,NFFT/2+1);
subplot(3,2,2)
plot(f,2*abs(Y(1:NFFT/2+1)))
xlabel('(Hz)')
ylabel('')
grid on
fn=1002;
fp=[80,230];
fs=[60,300];
Rp=3;
Rs=30;
Wp=fp/(fn/2);
Ws=fs/(fn/2);
[n,Wn]=buttord(Wp,Ws,Rp,Rs);
[b,a]=butter(n,Wn);
[H,F]=freqz(b,a,501,1002);
%plot(F,20*log10(abs(H)))
%axis([0 400 -30 3])
%xlabel('???? (Hz)'); ylabel('・ù??(dB)')
%grid on
%subplot(3,2,3)
pha=angle(H)*180/pi;
% plot(F,pha)
% xlabel('???? (Hz)'); ylabel('?à??')
% grid on
%±???????????
A=2*abs(Y(1:NFFT/2+1)).*(abs(H)');
subplot(3,2,3)
plot(f,A)
grid on
%×?ifft
for i=1:L-(NFFT/2+1)
A(1,((NFFT/2+1)+i))=A(1,((NFFT/2+1)-i));
end
z=A.*exp(1i*angle(Y));
z2=real(ifft(z));
subplot(3,2,4)
y=0.7*cos(2*pi*100*t)+sin(2*pi*200*t);
plot(t,y);
% subplot(3,2,5)
% plot(t,z2*NFFT/2);
%high filter---------------------------------
fn=1002;
fp=[60,300];
fs=[70,250];
Rp=3;
Rs=30;
Ws=fs/(fn/2);
[n,Wn]=buttord(Wp,Ws,Rp,Rs);
[b,a]=butter(n,Wn,'stop');
[H,F]=freqz(b,a,501,1002);
subplot(3,2,5)
plot(F,20*log10(abs(H)))
axis([0 400 -35 3])
xlabel('???? (Hz)'); ylabel('・ù??(dB)')
grid on
%subplot(3,2,6)
pha=angle(H)*180/pi;
plot(F,pha)
xlabel('???? (Hz)'); ylabel('?à??')
grid on
%-------------------
A=(2*abs(Y(1:NFFT/2+1))).*(abs(H)');
%ifft-----------------------------------
for i=1:L-(NFFT/2+1)
A(1,((NFFT/2+1)+i))=A(1,((NFFT/2+1)-i));
end
z=A.*exp(1i*angle(Y));
z2=real(ifft(z));
%subplot(3,2,6)
y=0.4*sin(2*pi*40*t)+5*cos(2*pi*400*t);
plot(t,y);
subplot(2,1,2)
plot(t,z2*NFFT/2);
