数字信号处理大作业代码及报告.zip

clear;close all;clc; %====原始语音信号=请自行设置文件地址audioread.m第74行filename = ('D:\pei.wav');% figure(1); filename = 'pei.wav'; [music,Fs] = audioread('pei.wav'); %模数转换 对语音信号进行采样 y1=music(:,1); %取单声道 y2=music(:,2); % sound(y1,Fs); Ns=length(y1); %采样点数Ns t=(0:Ns-1)/Fs; subplot(231); plot(t,y1);title('原始音频信号波形');xlabel('time/s');ylabel('幅度'); Y1=fft(y1); %对y1做快速傅里叶变换 magY1=abs(Y1); %取绝对值 即模值 magY1(1)=magY1(1)/Ns; %0频率点的信号幅值=fft变换后的幅值/Ns magY1(2:Ns)=magY1(2:Ns)/Ns*2; %其他频率点信号的幅值=fft变换后的幅值*2/Ns subplot(234); f=(0:Ns-1)*Fs/Ns; %Fs/Ns为频率分辨率 fft的每一点乘以分辨率代表此点的频率 plot(f,magY1);title('幅频特性曲线');xlabel('frequency/Hz');ylabel('幅度');axis([0 1e+4 0 0.03]); %================噪声====================================% t=(0:Ns-1)/Fs; noise=sin(2*pi*7500*t); %噪声为正弦波 幅值为1 频率为7500Hz subplot(232); plot(t,noise);title('噪声信号波形');xlabel('time/s');ylabel('幅度'); Noise=fft(noise); magN=abs(Noise); magN(1)=magN(1)/Ns; magN(2:Ns)=magN(2:Ns)/Ns*2; subplot(235); f=(0:length(Noise)-1)*Fs/Ns; plot(f,magN);title('幅频特性曲线');xlabel('frequency/Hz');ylabel('幅度');axis([0 1e+4 0 0.03]); %===============加性噪声================================% y1n=y1+noise'; % sound(y1n,Fs); t=(0:length(y1n)-1)/Fs; subplot(233); plot(t,y1n);title('加噪后音频信号波形');xlabel('time/s');ylabel('幅度'); Xn=fft(y1n); magXn=abs(Xn); magXn(1)=magXn(1)/Ns; magXn(2:Ns)=magXn(2:Ns)*2/Ns; subplot(236); f=(0:length(Xn)-1)*Fs/length(Xn); plot(f,magXn);title('幅频特性曲线');xlabel('frequency/Hz');ylabel('幅度');axis([0 1e+4 0 0.03]); %==================巴特沃斯低通滤波器======================% fp=4900;fs=7450;As=100;Rp=1; %缺点：过渡带要求较宽 阶数高 26阶 OmegaP=2*pi*fp;OmegaS=2*pi*fs; %模拟角频率 wp=OmegaP/Fs;ws=OmegaS/Fs; %数字频率 wp1=2*Fs*tan(wp/2); %频率预矫正 ws1=2*Fs*tan(ws/2); [N,wc]=buttord(wp1,ws1,Rp,As,'s'); %切比雪夫II型阶数与截止频率 [b,a]=butter(N,wc,'s'); %模拟切比雪夫滤波器分子、分母系数 [bz,az]=bilinear(b,a,Fs); %双线性变换AF到DF w0=[wp,ws]; Hx1=freqz(bz,az,w0); %计算在矢量w0中指定频率处的得频率响应 [H,w]=freqz(bz,az); dbHx1=-20*log10(abs(Hx1)/max(abs(H))); %求滤波器截止频率处的衰减 figure(4); plot(w/(2*pi)*Fs/1000,20*log10(abs(H))); title('巴特沃斯低通数字滤波器幅频特性'); xlabel('f(KHz)');ylabel('dB'); %----滤波后音频信号----% figure(2); out=filter(bz,az,y1n); t=(0:length(out)-1)/Fs; subplot(221); plot(t,out);title('巴特沃斯低通滤波后音频信号波形');xlabel('time/s');ylabel('幅度'); OUT=fft(out); magOUT=abs(OUT); magOUT(1)=magOUT(1)/Ns; magOUT(2:Ns)=magOUT(2:Ns)/Ns*2; subplot(223); f=(0:length(OUT)-1)*Fs/length(OUT); plot(f,magOUT);title('幅频特性曲线');xlabel('frequency/Hz');ylabel('幅度');axis([0 1e+4 0 0.03]); sound(out,Fs); %==================切比雪夫II型低通滤波器======================% fp=6900;fs=7450;As=100;Rp=1; %通带截止频率与阻带截止频率的确定采用逼近的方法，在满足阻带最小衰减As前提下逐渐缩小过渡带 OmegaP=2*pi*fp;OmegaS=2*pi*fs; %模拟角频率 wp=OmegaP/Fs;ws=OmegaS/Fs; %数字频率 wp1=2*Fs*tan(wp/2); %频率预矫正 ws1=2*Fs*tan(ws/2); [N,wc]=cheb2ord(wp1,ws1,Rp,As,'s'); %切比雪夫II型阶数与截止频率 30阶 [b,a]=cheby2(N,As,wc,'s'); %模拟切比雪夫滤波器分子、分母系数 [bz,az]=bilinear(b,a,Fs); %双线性变换AF到DF figure(5); w0=[wp,ws]; Hx2=freqz(bz,az,w0); %计算在矢量w0中指定频率处的得频率响应 [H,w]=freqz(bz,az); dbHx2=-20*log10(abs(Hx2)/max(abs(H))); %求滤波器截止频率处的衰减 figure(5); plot(w/(2*pi)*Fs/1000,20*log10(abs(H))); title('切比雪夫II型低通数字滤波器幅频特性'); xlabel('f(KHz)');ylabel('dB'); %-----------------滤波后音频信号-------------% figure(2); out=filter(bz,az,y1n); t=(0:length(out)-1)/Fs; subplot(222);plot(t,out);title('切比雪夫II型低通滤波后音频信号波形');xlabel('time/s');ylabel('幅度'); OUT=fft(out); magOUT=abs(OUT); magOUT(1)=magOUT(1)/Ns; magOUT(2:Ns)=magOUT(2:Ns)/Ns*2; subplot(224); f=(0:length(OUT)-1)*Fs/length(OUT); plot(f,magOUT);title('幅频特性曲线');xlabel('frequency/Hz');ylabel('幅度');axis([0 1e+4 0 0.03]); % sound(out,Fs); %==================巴特沃斯带阻滤波器======================% fp1=7000;fs1=7450;fs2=7550;fp2=8000;As=100;Rp=1; %通带截止频率与阻带截止频率的确定采用逼近的方法，在满足阻带最小衰减As前提下逐渐缩小过渡带 OmegaP1=2*pi*fp1;OmegaS1=2*pi*fs1; %模拟角频率 OmegaP2=2*pi*fp2;OmegaS2=2*pi*fs2; wp1=OmegaP1/Fs;ws1=OmegaS1/Fs; %数字频率 wp2=OmegaP2/Fs;ws2=OmegaS2/Fs; wp11=2*Fs*tan(wp1/2); %频率预畸 ws11=2*Fs*tan(ws1/2); wp22=2*Fs*tan(wp2/2); ws22=2*Fs*tan(ws2/2); wp=[wp11,wp22]; ws=[ws11,ws22]; [N,wc]=buttord(wp,ws,Rp,As,'s'); %确定滤波器阶数与截止频率 6阶 [b,a]=butter(N,wc,'stop','s'); %模拟滤波器分子、分母系数 [bz,az]=bilinear(b,a,Fs); %双线性变换AF到DF w0=[wp1,ws1,ws2,wp2]; Hx3=freqz(bz,az,w0); %计算在矢量w0中指定频率处的得频率响应 [H,w]=freqz(bz,az); dbHx3=-20*log10(abs(Hx3)/max(abs(H))); %求滤波器截止频率处的衰减 figure(6); plot(w/(2*pi)*Fs/1000,20*log10(abs(H))); title('巴特沃斯带阻数字滤波器幅频特性'); xlabel('f(KHz)');ylabel('dB'); %----滤波后音频信号----% figure(3); out1=filter(bz,az,y1n); t=(0:length(out1)-1)/Fs; subplot(221); plot(t,out1);title('巴特沃斯带阻滤波后音频信号波形');xlabel('time/s');ylabel('幅度'); OUT1=fft(out1); magOUT1=abs(OUT1); magOUT1(1)=magOUT1(1)/Ns; magOUT1(2:Ns)=magOUT1(2:Ns)/Ns*2; subplot(223); f=(0:length(OUT1)-1)*Fs/length(OUT1); plot(f,magOUT1);title('幅频特性曲线');xlabel('frequency/Hz');ylabel('幅度');axis([0 1e+4 0 0.03]); % sound(out1,Fs); %==================切比雪夫II型带阻滤波器======================% fp1=6800;fs1=7450;fs2=7550;fp2=7800;As=100;Rp=1; OmegaP1=2*pi*fp1;OmegaS1=2*pi*fs1; %模拟角频率 OmegaP2=2*pi*fp2;OmegaS2=2*pi*fs2; wp1=OmegaP1/Fs;ws1=OmegaS1/Fs; %数字频率 wp2=OmegaP2/Fs;ws2=OmegaS2/Fs; wp11=2*Fs*tan(wp1/2); %频率预畸 ws11=2*Fs*tan(ws1/2); wp22=2*Fs*tan(wp2/2); ws22=2*Fs*tan(ws2/2); wp=[wp11,wp22]; ws=[ws11,ws22]; [N,wc]=cheb2ord(wp,ws,Rp,As,'s'); %切比雪夫II型阶数与截止频率 6阶 [b,a]=cheby2(N,As,wc,'stop','s'); %模拟切比雪夫滤波器分子、分母系数 [bz,az]=bilinear(b,a,Fs); %双线性变换AF到DF w0=[wp1,ws1,ws2,wp2]; Hx4=freqz(bz,az,w0); %计算在矢量w0中指定频率处的得频率响应 [H,w]=freqz(bz,az); dbHx4=-20*log10(abs(Hx4)/max(abs(H))); %求滤波器截止频率处的衰减 figure(7); plot(w/(2*pi)*Fs/1000,20*log10(abs(H))); title('切比雪夫II型带阻数字滤波器幅频特性'); xlabel('f(KHz)');ylabel('dB'); %----滤波后音频信号----% figure(3); out1=filter(bz,az,y1n); t=(0:length(out1)-1)/Fs; subplot(222);plot(t,out1);title('切比雪夫II型带阻滤波后音频信号波形');xlabel('time/s');ylabel(