用matlab设计滤波器

这个信号的频率分量分别为30、150和600Hz，因此可分别设计一个低通、带通和高通的滤波器来提取。以FIR滤波器为例，程序如下：

clearfs=2000t=(1:1000)/fs

x=10*cos(2*pi*30*t)+cos(2*pi*150*t)+5*cos(2*pi*600*t)

L=length(x)N=2^(nextpow2(L))Hw=fft(x,N)

figure(1)subplot(2,1,1)plot(t,x)

grid ontitle('滤波前信号x')xlabel('时间/s')% 原始信号

subplot(2,1,2)plot((0:N-1)*fs/L,abs(Hw))% 查看信号频谱

grid ontitle('滤波前信号频谱图')xlabel('频率/Hz')ylabel('振幅|H(e^jw)|')

%% x_1=10*cos(2*pi*30*t)

Ap=1As=60% 定义通带及阻带衰减

dev=[(10^(Ap/20)-1)/(10^(Ap/20)+1),10^(-As/20)]% 计算偏移量

mags=[1,0]% 低通

fcuts=[60,100]% 边界频率

[N,Wn,beta,ftype]=kaiserord(fcuts,mags,dev,fs)% 估算FIR滤波器阶数

hh1=fir1(N,Wn,ftype,kaiser(N+1,beta))% FIR滤波器设计

x_1=filter(hh1,1,x)% 滤波

x_1(1:ceil(N/2))=[]% 群延时N/2，删除无用信号部分

L=length(x_1)N=2^(nextpow2(L))Hw_1=fft(x_1,N)

figure(2)subplot(2,1,1)plot(t(1:L),x_1)

grid ontitle('x_1=10*cos(2*pi*30*t)')xlabel('时间/s')

subplot(2,1,2)plot((0:N-1)*fs/L,abs(Hw_1))% 查看信号频谱

grid ontitle('滤波后信号x_1频谱图')xlabel('频率/Hz')ylabel('振幅|H(e^jw)|')

%% x_2=cos(2*pi*150*t)

Ap=1As=60% 定义通带及阻带衰减

dev=[10^(-As/20),(10^(Ap/20)-1)/(10^(Ap/20)+1),10^(-As/20)]% 计算偏移量

mags=[0,1,0]% 带通

fcuts=[80,120,180,220]% 边界频率

[N,Wn,beta,ftype]=kaiserord(fcuts,mags,dev,fs)% 估算FIR滤波器阶数

hh2=fir1(N,Wn,ftype,kaiser(N+1,beta))% FIR滤波器设计

x_2=filter(hh2,1,x)% 滤波

x_2(1:ceil(N/2))=[]% 群延时N/2，删除无用信号部分

L=length(x_2)N=2^(nextpow2(L))Hw_2=fft(x_2,N)

figure(3)subplot(2,1,1)plot(t(1:L),x_2)

grid ontitle('x_2=cos(2*pi*150*t)')xlabel('时间/s')

subplot(2,1,2)plot((0:N-1)*fs/L,abs(Hw_2))% 查看信号频谱

grid ontitle('滤波后信号x_2频谱图')xlabel('频率/Hz')ylabel('振幅|H(e^jw)|')

%% x_3=5*cos(2*pi*600*t)

Ap=1As=60% 定义通带及阻带衰减

dev=[10^(-As/20),(10^(Ap/20)-1)/(10^(Ap/20)+1)]% 计算偏移量

mags=[0,1]% 高通

fcuts=[500,550]% 边界频率

[N,Wn,beta,ftype]=kaiserord(fcuts,mags,dev,fs)% 估算FIR滤波器阶数

hh2=fir1(N,Wn,ftype,kaiser(N+1,beta))% FIR滤波器设计

x_3=filter(hh2,1,x)% 滤波

x_3(1:ceil(N/2))=[]% 群延时N/2，删除无用信号部分

L=length(x_3)N=2^(nextpow2(L))Hw_3=fft(x_3,N)

figure(4)subplot(2,1,1)plot(t(1:L),x_3)

grid ontitle('x_3=5*cos(2*pi*600*t)')xlabel('时间/s')

subplot(2,1,2)plot((0:N-1)*fs/L,abs(Hw_3))% 查看信号频谱

grid ontitle('滤波后信号x_3频谱图')xlabel('频率/Hz')ylabel('振幅|H(e^jw)|')

你自己整合吧，我没时间帮你整合，我给你提供一些程序：

绝对正确的代码：程序1：

fs=22050 %语音信号采样频率为22050

x1=wavread('Windows Critical Stop.wav')%读取语音信号的数据，赋给变量x1

sound(x1,22050) %播放语音信号

y1=fft(x1,1024) %对信号做1024点FFT变换

f=fs*(0:511)/1024

figure(1)

plot(x1) %做原始语音信号的时域图形

title('原始语音信号')

xlabel('time n')

ylabel('fuzhi n')

figure(2)

freqz(x1) %绘制原始语音信号的频率响应图

title('频率响应图')

figure(3)

subplot(2,1,1)

plot(abs(y1(1:512))) %做原始语音信号的FFT频谱图

title('原始语音信号FFT频谱')

subplot(2,1,2)

plot(f,abs(y1(1:512)))

title('原始语音信号频谱')

xlabel('Hz')

ylabel('fuzhi')

程序2：

fs=22050 %语音信号采样频率为22050

x1=wavread('Windows Critical Stop.wav')%读取语音信号的数据，赋给变量x1

t=0:1/22050:(size(x1)-1)/22050

y1=fft(x1,1024) %对信号做1024点FFT变换

f=fs*(0:511)/1024

x2=randn(1,length(x1)) %产生一与x长度一致的随机信号

sound(x2,22050)

figure(1)

plot(x2) %做原始语音信号的时域图形

title('高斯随机噪声')

xlabel('time n')

ylabel('fuzhi n')

randn('state',0)

m=randn(size(x1))

x2=0.1*m+x1

sound(x2,22050)%播放加噪声后的语音信号

y2=fft(x2,1024)

figure(2)

plot(t,x2)

title('加噪后的语音信号')

xlabel('time n')

ylabel('fuzhi n')

figure(3)

subplot(2,1,1)

plot(f,abs(y2(1:512)))

title('原始语音信号频谱')

xlabel('Hz')

ylabel('fuzhi')

subplot(2,1,2)

plot(f,abs(y2(1:512)))

title('加噪后的语音信号频谱')

xlabel('Hz')

ylabel('fuzhi')

根据以上代码，你可以修改下面有错误的代码

程序3:双线性变换法设计Butterworth滤波器

fs=22050

x1=wavread('h:\课程设计2\shuzi.wav')

t=0:1/22050:(size(x1)-1)/22050

Au=0.03

d=[Au*cos(2*pi*5000*t)]'

x2=x1+d

wp=0.25*pi

ws=0.3*pi

Rp=1

Rs=15

Fs=22050

Ts=1/Fs

wp1=2/Ts*tan(wp/2)%将模拟指标转换成数字指标

ws1=2/Ts*tan(ws/2)

[N,Wn]=buttord(wp1,ws1,Rp,Rs,'s') %选择滤波器的最小阶数

[Z,P,K]=buttap(N) %创建butterworth模拟滤波器

[Bap,Aap]=zp2tf(Z,P,K)

[b,a]=lp2lp(Bap,Aap,Wn)

[bz,az]=bilinear(b,a,Fs) %用双线性变换法实现模拟滤波器到数字滤波器的转换

[H,W]=freqz(bz,az)%绘制频率响应曲线

figure(1)

plot(W*Fs/(2*pi),abs(H))

grid

xlabel('频率／Hz')

ylabel('频率响应幅度')

title('Butterworth')

f1=filter(bz,az,x2)

figure(2)

subplot(2,1,1)

plot(t,x2) %画出滤波前的时域图

title('滤波前的时域波形')

subplot(2,1,2)

plot(t,f1)%画出滤波后的时域图

title('滤波后的时域波形')

sound(f1,22050) %播放滤波后的信号

F0=fft(f1,1024)

f=fs*(0:511)/1024

figure(3)

y2=fft(x2,1024)

subplot(2,1,1)

plot(f,abs(y2(1:512)))%画出滤波前的频谱图

title('滤波前的频谱')

xlabel('Hz')

ylabel('fuzhi')

subplot(2,1,2)

F1=plot(f,abs(F0(1:512))) %画出滤波后的频谱图

title('滤波后的频谱')

xlabel('Hz')

ylabel('fuzhi')

程序4:窗函数法设计滤波器：

fs=22050

x1=wavread('h:\课程设计2\shuzi.wav')

t=0:1/22050:(size(x1)-1)/22050

Au=0.03

d=[Au*cos(2*pi*5000*t)]'

x2=x1+d

wp=0.25*pi

ws=0.3*pi

wdelta=ws-wp

N=ceil(6.6*pi/wdelta) %取整

wn=(0.2+0.3)*pi/2

b=fir1(N,wn/pi,hamming(N+1)) %选择窗函数，并归一化截止频率

figure(1)

freqz(b,1,512)

f2=filter(bz,az,x2)

figure(2)

subplot(2,1,1)

plot(t,x2)

title('滤波前的时域波形')

subplot(2,1,2)

plot(t,f2)

title('滤波后的时域波形')

sound(f2,22050) %播放滤波后的语音信号

F0=fft(f2,1024)

f=fs*(0:511)/1024

figure(3)

y2=fft(x2,1024)

subplot(2,1,1)

plot(f,abs(y2(1:512)))

title('滤波前的频谱')

xlabel('Hz')

ylabel('fuzhi')

subplot(2,1,2)

F2=plot(f,abs(F0(1:512)))

title('滤波后的频谱')

xlabel('Hz')

ylabel('fuzhi')

将模拟频率转化为数字频率，设取样时间为T（要满足抽样定理） Ωp=2π*fp*T Ωs=2π*fs*T 过渡带宽度△Ω=Ωp-Ωs 阻带衰减已经超过74db，要选用Kaiser窗了，Kaiser的参数可变，要根据公式确定滤波器的参数一般都选用Ⅰ型线性相位滤波器即滤波器阶数M为偶数，程序如下： wp=ws=Ap=1As=100Rp=1-10.^(-0.05*Ap)Rs=10.^(-0.05*As)f=[fp fs]a=[0 1]dev=[Rp Rs][M,wc,beta,ftype]=kaiserord(f,a,dev)M=mod(M,2)+Mh=fir1(M,wc,ftype,kaiser(M+1,beta))omega=linspace(0,pi,512)mag=freqz(h,[1],omega)plot(omega/pi,20*log10(abs(mag)))gridomega1=linspace(0,wp,512)h1=freqz(h,[1],omega1)omega2=linspace(ws,pi,512)h2=freqz(h,[1],omega2)fprintf('Ap=%.4f\n',-20*log10(min(abs(h1))))fprintf('As=%.4f\n',-20*log10(max(abs(h2))))运行程序可以得到滤波器的通阻带衰减，画出频率响应，若同阻带衰减不满足要求还可以使用滤波器的优化，一般使用的等波纹FIR进行优化

---