数字滤波器与音乐信号滤波

数字滤波器与音乐信号滤波报告1数字滤波器与音乐信号滤波【目的】(1)掌握IIR和FIR数字滤波器的设计指标、设计方法及各自的特点。(2)掌握BW型、CBI型和椭圆型IIR数字滤波器的特点。(3)掌握窗函数对FIR数字滤波器设计的影响。(4)培养学生自主学习能力，以及发现问题、分析问题和解决问题的能力。【知识点】数字滤波器指标IIR数字滤波器FIR数字滤波器【背景知识】音乐中的和弦一般由3个（也可大于3个）具有一定规律的乐音构成，D大调和弦由DbGA三个乐音构成，即由频率为293.66,369.99,440Hz的正弦信号组成。【研讨内容】Dhexian.wav是将D大调和弦以8000Hz抽样所得的数字音乐信号，试设计数字滤波器从和弦中分离出乐音A，即440Hz的频率分量。分离乐音A的模拟滤波器的设计指标如下：模拟滤波器指标1：fp=440Hz，Ap=1dB，fs=370Hz，As=20dB；模拟滤波器指标2：fp=440Hz，Ap=1dB，fs=370Hz，As=40dB；模拟滤波器指标3：fp=440Hz，Ap=1dB，fs=370Hz，As=50dB；模拟滤波器指标4：fp=440Hz，Ap=1dB，fs=380Hz，As=50dB；模拟滤波器指标5：fp=440Hz，Ap=1dB，fs=360Hz，As=50dB。要求：1.设计IIR数字滤波器（1）确定能够取代题目中5个模拟滤波器的IIR数字滤波器的设计指标。数字滤波器指标1：Ωp=0.11πrad，Ap=1dB，Ωs=0.0925πrad，As=20dB；数字滤波器指标2：Ωp=0.11πrad，Ap=1dB，Ωs=0.0925πrad，As=40dB；数字滤波器指标3：Ωp=0.11πrad，Ap=1dB，Ωs=0.0925πrad，As=50dB；数字滤波器指标4：Ωp=0.11πrad，Ap=1dB，Ωs=0.095πrad，As=50dB；数字滤波器指标5：Ωp=0.11πrad，Ap=1dB，Ωs=0.09πrad，As=50dB。（2）设计满足上述指标的CBI型IIR数字滤波器，综合考虑信号处理成本和本组成员对滤波效果的认可，你们组认为和弦中分离出乐音A的数字滤波器指标选择哪个比较合适？编写matlab程序如下：clear;closeall;fp=440;fs=370;Ap=1;As=20;%模拟滤波器指标1%fp=440;fs=370;Ap=1;As=40;%模拟滤波器指标2%fp=440;fs=370;Ap=1;As=50;%模拟滤波器指标3%fp=440;fs=380;Ap=1;As=50;%模拟滤波器指标4%fp=440;fs=360;Ap=1;As=50;%模拟滤波器指标5Wp=(fp*2*pi)/8000;Ws=(fs*2*pi)/8000;数字滤波器与音乐信号滤波报告2wp=Wp/pi;ws=Ws/pi;%CBI型IIR数字滤波器[N,wc]=cheb1ord(wp,ws,Ap,As);[num,den]=cheby1(N,Ap,wc,'high');%figure%zplane(num,den);%零极点分布图w=linspace(0,pi,1024);h=freqz(num,den,w);figure('name','滤波器增益');plot(w/pi,20*log10(abs(h)));gridon;%滤波器增益figure('name','滤波器幅频特性与相频特性');subplot(2,1,1);plot(w/pi,abs(h));gridon;%滤波器幅度响应title('幅度响应');ylabel('幅度谱');xlabel('归一化角频率');subplot(2,1,2);plot(w/pi,angle(h));gridon;%滤波器相位响应title('相位响应');ylabel('相位谱');xlabel('归一化角频率');w=[wpws]*pi;h=freqz(num,den,w);fprintf('Ap=%.4f\n',-20*log10(abs(h(1))));fprintf('As=%.4f\n',-20*log10(abs(h(2))));%读取音乐获得离散数据点和采样频率[data,fs]=audioread(Ddiaohexian.wav);L=length(data);%获得时域抽样点数Ny1=data(:,1);y2=filter(num,den,y1);sound(y2,fs);X1=fftshift(fft(y1,L));X2=fftshift(fft(y2,L));figure('name','原时域波形');plot(y1);%原时域波形figure('name','滤波时域波形');plot(y2);%滤波时域波形f=(-fs/2+(0:L-1)*fs/L);figure('name','原幅度频谱');plot(f,abs(X1));gridon;%原幅度频谱title('原信号幅度谱');ylabel('幅度谱');xlabel('频率(Hz)');figure('name','滤波幅度频谱');plot(f,abs(X2));gridon;%滤波幅度频谱title('滤波后幅度谱');数字滤波器与音乐信号滤波报告3ylabel('幅度谱');xlabel('频率(Hz)');原信号幅度频谱：设计指标1：滤波器阶数6，As=20.0157dB滤波器幅度响应和相位响应：数字滤波器与音乐信号滤波报告4信号滤波后频谱：设计指标2：滤波器阶数10，As=41.2039dB滤波器幅度响应和相位响应：信号滤波后频谱：数字滤波器与音乐信号滤波报告5设计指标3：滤波器阶数12，As=51.8221dB滤波器幅度响应和相位响应：信号滤波后频谱：数字滤波器与音乐信号滤波报告6设计指标4：滤波器阶数13，As=51.3296dB滤波器幅度响应和相位响应：信号滤波后频谱：数字滤波器与音乐信号滤波报告7设计指标5：滤波器阶数11，As=51.2348dB滤波器幅度响应和相位响应：信号滤波后频谱：数字滤波器与音乐信号滤波报告8总结由各指标设计的滤波器如下：设计指标CBI型IIR数字滤波器阶数As1620.015721041.203931251.822141351.329651151.2348分析：考虑成本及滤波效果，选择设计指标2的CBI型IIR数字滤波器分离乐音A。（3）设计能够取代满足指标1的模拟滤波器的BW型、CBI型、CBII型和椭圆型IIR数字滤波器，并画出其零极点分布图。编写matlab程序如下：clear;closeall;fp=440;fs=370;Ap=1;As=20;Wp=(fp*2*pi)/8000;Ws=(fs*2*pi)/8000;wp=Wp/pi;ws=Ws/pi;%BW型滤波器[N,wc]=buttord(wp,ws,Ap,As);[num,den]=butter(N,wc,'high');figurezplane(num,den);%零极点分布图w=linspace(0,pi,1024);数字滤波器与音乐信号滤波报告9h=freqz(num,den,w);figure('name','滤波器增益');plot(w/pi,20*log10(abs(h)));gridon;%滤波器增益figure('name','滤波器幅频特性与相频特性');subplot(2,1,1);plot(w/pi,abs(h));gridon;%滤波器幅度响应title('幅度响应');ylabel('幅度谱');xlabel('归一化角频率');subplot(2,1,2);plot(w/pi,angle(h));gridon;%滤波器相位响应title('相位响应');ylabel('相位谱');xlabel('归一化角频率');w=[wpws]*pi;h=freqz(num,den,w);fprintf('Ap=%.4f\n',-20*log10(abs(h(1))));fprintf('As=%.4f\n',-20*log10(abs(h(2))));%读取音乐获得离散数据点和采样频率[data,fs]=audioread(Ddiaohexian.wav);L=length(data);%获得时域抽样点数Ny1=data(:,1);y2=filter(num,den,y1);sound(y2,fs);X1=fftshift(fft(y1,L));X2=fftshift(fft(y2,L));figure('name','原时域波形');plot(y1);%原时域波形figure('name','滤波时域波形');plot(y2);%滤波时域波形f=(-fs/2+(0:L-1)*fs/L);figure('name','原幅度频谱');plot(f,abs(X1));gridon;%原幅度频谱title('原信号幅度谱');ylabel('幅度谱');xlabel('频率(Hz)');figure('name','滤波幅度频谱');plot(f,abs(X2));gridon;%滤波幅度频谱title('滤波后幅度谱');ylabel('幅度谱');xlabel('频率(Hz)');BW型：滤波器阶数17，Ap=0.9616dB，As=20.0025dB滤波器幅度响应和相位响应：数字滤波器与音乐信号滤波报告10信号滤波后频谱：零极点分布图：数字滤波器与音乐信号滤波报告11CBI型：1234561234560.38512.31065.77657.70205.77652.31060.3851()14.26837.85827.96554.72851.58730.2458zzz+zz+zHzzzzzz+z−−−−−−−−−−−−−+−−=−+−+−滤波器阶数6，Ap=1.0000dB，As=20.0157dB滤波器幅度响应和相位响应、信号滤波后频谱同（2）。零极点分布图：数字滤波器与音乐信号滤波报告12CBII型：1234561234560.74714.388110.831514.380910.83154.38810.7471()15.299611.838714.25889.76043.59860.5581zzz+zz+zHzzzzzz+z−−−−−−−−−−−−−+−−=−+−+−滤波器阶数6，Ap=0.9967dB，As=20.0000dB滤波器幅度响应和相位响应：信号滤波后频谱：数字滤波器与音乐信号滤波报告13零极点分布图：椭圆型：123412340.67932.63363.90.99730.3582()13.31874.31222.58290.6080zzzzHzzzz+z−−−−−−−−−+−+=−+−滤波器阶数4，Ap=1.0000dB，As=22.0230dB滤波器幅度响应和相位响应：数字滤波器与音乐信号滤波报告14信号滤波后频谱：零极点分布图：（4）比较(3)中不同滤波器的阶数以及零极点分布特点。阶数：椭圆型4阶CBII型6阶=CBI型6阶BW型17阶，零极点分布全部在右半平面，且关于实轴对称。数字滤波器与音乐信号滤波报告152.设计FIR数字滤波器（1）利用窗函数法设计FIR数字滤波器，确定能够取代题目中5个模拟滤波器的数字滤波器参数，即待逼近理想滤波器的截止频率、窗函数的类型、滤波器的阶数。编写matlab程序如下：clear;closeall;fp=440;fs=370;Ap=1;As=20;%fp=440;fs=370;Ap=1;As=40;%fp=440;fs=370;Ap=1;As=50;%fp=440;fs=380;Ap=1;As=50;%fp=440;fs=360;Ap=1;As=50;wp=(fp*2*pi)/8000;ws=(fs*2*pi)/8000;%读取音乐获得离散数据点和采样频率[data,fs]=audioread(Ddiaohexian.wav)