MATLAB 数字信号处理 滤波器设计

[帅小龙]

设计目的和内容 通过数字信号处理的大作业设计，巩固和运用数字信号处理课程中的理论知识和实验技能，掌握最基本的数字信号处理的理论和方法。本次设计的题目是：语音信号的处理与滤波。 具体来说就是对一段语音信号进行采样；画出采样后语音信号的时域波形和频谱图；然后给采样信号加以噪声，采用双线性变换设计IIR滤波器；然后用设计的滤波器对加噪声的信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；分别回放原始、加噪声、滤波后语音信号，感受其直观变化。 设计原理 1.语音信号的采集 语音信号是一种模拟信号，首先须经过采样将其转换为数字信号，实质是把连续信号变为脉冲或数字序列。我们用Windows10自带的录音机录制一段mp3格式的音频。然后用MATLAB的audioread函数采样。配合sound函数来使用。 2.语音信号的频谱特性 在MATLAB中，可以利用函数fft对信号进行快速傅里叶变换，得到信号的频谱特性。 3.语音信号的时域分析 语音信号的时域分析就是分析和提取语音信号的时域参数。。进行语音分析时，最先接触到并且也是最直观的是它的时域波形。语音信号本身就是时域信号，因而时域分析是最早使用，也是应用最广泛的一种分析方法，这种方法直接利用语音信号的时域波形。 设计过程 1.语音信号的采集 我们用Windows10自带的录音机录制一段mp3格式的音频，我录制了一段30秒左右的音乐。将然后用MATLAB的audioread函数对这一段语音进行采样。配合sound函数使用可播放该语音信号 2.对语音信号进行频谱分析 处理语音信号的时域波形图，对语音信号进行快速傅立叶变换，得到信号的频谱。 运行后会播放原始音乐，并且出现的时域波形图和频谱图： 语音信号加噪声并回放加噪声语音信号 使用randn函数产生噪声，并将其加入到原始语音信号中。 运行后，会播放加噪声后的语音信号，可以听见有明显的“雪花”声。并且出现如下的时域波形图和频谱图： 根据指标设计数字滤波器 可以运用窗函数法、双线性变换法设计低通、高通、带通三种滤波器。我选用双线性变换法进行设计低通滤波器。 运行后可以得到滤波器的特性图 用滤波器对信号进行滤波并回放滤波语音信号 运行后，噪声比加噪音的信号要小，但是还能听见一些“雪花”声。其波形，频谱如下： （我还使用了ellipord函数，分别设计了低通，高通，带通滤波器，但是除了低通滤波器具有一定的效果以外，高通和带通滤波器都没有很好地效果） 设计总结 通过这次课程设计，让我对MATLAB有了深入的了解。由于电脑只安装了2018a版本的MATLAB。于我上课所用的MATLAB具有一定不同，所以去查了许多的资料，才完成了这次的大作业。这次设计主要参考了CSDN上的一些大佬做的设计，看了他们很多的代码，自己才慢慢的懂得了如何去写一些代码，和一些代码的作用。这次设计让我懂得了，要完成一件事需要很多很多的努力。努力就会有收获。 附录：实验程序 双击图标打开录音，程序运行中按任意键继续 ===============大作业实验程序=============== [x,fs]=audioread('E:\MATLAB\music.mp3');%声音读取 sound(x,fs); %播放原始的语音 x=x(:,1); FS=length(x);%求出原始语音信号长度 X=fft(x);%傅里叶变换 t=(0:FS-1)/fs; figure(1) subplot(211); plot(t,x); title('原始语音信号时域波形'); xlabel('时间'); ylabel('幅度'); grid on;%显示坐标轴网格线 subplot(212); plot(abs(X)); title('原始语音信号频谱'); xlabel('频率'); ylabel('幅度'); axis([0 1400000 0 20]); grid on; pause;%暂停 FS=length(x); noise=0.01*randn(FS,1);%加噪声 s=x+noise; sound(s,fs); %播放加噪的语音 N=length(s); %画出加噪之后，其时域频域 S=fft(s,N); figure(2) subplot(211) plot(s); title('加噪声后信号波形') xlabel('时间'); ylabel('幅度'); grid on; subplot(212) plot(abs(fftshift(S))); xlabel('频率'); ylabel('幅度'); axis([0 1400000 0 20]); title('加噪声后信号信号频谱'); grid on; pause; %设计IIR低通滤波器 Rp=2;Rs=80; Ft=8000;%采样频率 Fp=2000; Fs=1000; wp=2*pi*Fp/Ft; ws=2*pi*Fs/Ft;%求出待设计的模拟滤波器的边界频率 [N,wn]=buttord(wp,ws,Rp,Rs,'s'); %设计过渡模拟滤波器:计算相应的模拟滤波器阶数N和截止频率wn [b,a]=butter(N,wn,'s'); %S域频率响应的参数即:计算相应的模拟滤波器系统函数 [bz,az]=bilinear(b,a,0.5); %利用双线性变换转换成数字滤波器 figure(3);%绘图 [h,w]=freqz(bz,az);%低通滤波器特性 title('IIR低通滤波器'); plot(w*fs/(2*pi),abs(h)); grid on %滤波 z=filter(bz,az,s); %滤波后画出其频谱，波形 Z=fft(z); %滤波后的信号频谱 figure(4) subplot(211); plot(z); title('低通滤波后的信号波形'); xlabel('时间'); ylabel('幅度'); grid on; subplot(212); plot(abs(fftshift(Z))); title('低通滤波后信号的频谱'); xlabel('频率'); ylabel('幅度'); axis([0 1400000 0 20]); grid on; sound(z,fs); %回放滤波后的信号 pause; =================额外的试验程序================ [x,fs]=audioread('E:\MATLAB\music.mp3');%声音读取 x=x(:,1); FS=length(x);%求出原始语音信号长度 X=fft(x);%傅里叶变换 t=(0:FS-1)/fs; figure(1) subplot(211); plot(t,x); title('原始语音信号时域波形'); xlabel('时间'); ylabel('幅度'); grid on;%显示坐标轴网格线 subplot(212); plot(abs(X)); title('原始语音信号频谱'); xlabel('频率'); ylabel('幅度'); axis([0 700000 0 20]); grid on; sound(x,fs); %播放原始的语音 pause;%暂停 FS=length(x); noise=0.01*randn(FS,1);%加噪声 s=x+noise; N=length(s); %画出加噪之后，其时域频域 S=fft(s,N); figure(2) subplot(211) plot(s); title('加噪声后信号波形') xlabel('时间'); ylabel('幅度'); grid on; subplot(212) plot(abs(fftshift(S))); xlabel('频率'); ylabel('幅度'); axis([0 700000 0 20]); title('加噪声后信号信号频谱'); grid on; sound(s,fs); %播放加噪的语音 pause; %设计IIR低通滤波器 Rp=2; Rs=80; Ft=8000;%采样频率 Fp=2000; Fs=1000; wp=2*pi*Fp/Ft; ws=2*pi*Fs/Ft;%求出待设计的模拟滤波器的边界频率 [N,wn]=buttord(wp,ws,Rp,Rs,'s'); [b,a]=butter(N,wn,'s'); [bz,az]=bilinear(b,a,0.5); figure(3);%绘图 [h,w]=freqz(bz,az);%低通滤波器特性 title('IIR低通滤波器'); plot(w*fs/(2*pi),abs(h)); grid on %滤波 z=filter(bz,az,s); %滤波后画出其频谱，波形 Z=fft(z); %滤波后的信号频谱 figure(4) subplot(211); plot(z); title('低通滤波后的信号波形'); xlabel('时间'); ylabel('幅度'); grid on; subplot(212); plot(abs(fftshift(Z))); title('低通滤波后信号的频谱'); xlabel('频率'); ylabel('幅度'); axis([0 700000 0 20]); grid on; sound(z,fs); %回放滤波后的信号 pause; %设计IIR高通滤波器 Rp=2; Rs=80; Ft=8000;%采样频率 Fp=1200; Fs=2200; wp=2*pi*Fp/Ft; ws=2*pi*Fs/Ft;%求出待设计的模拟滤波器的边界频率 [N,wn]=ellipord(wp,ws,Rp,Rs,'s'); [b,a]=ellip(N,Rp,Rs,wn,'high','s'); [bz,az]=bilinear(b,a,0.5); %绘图 figure(5); [h,w]=freqz(bz,az);%高通滤波器特性 title('IIR高通滤波器'); plot(w*fs/(2*pi),abs(h)); grid on %滤波 z=filter(bz,az,s); %滤波后画出其频谱，波形 Z=fft(z); %滤波后的信号频谱 figure(6) subplot(211); plot(z); title('高通滤波后的信号波形'); xlabel('时间'); ylabel('幅度'); grid on; subplot(212); plot(abs(fftshift(Z))); title('高通滤波后信号的频谱'); xlabel('频率'); ylabel('幅度'); axis([0 700000 0 20]); grid on; sound(z,fs); %回放滤波后的信号 pause; %设计IIR带通滤波器 Rp=2; Rs=80; Ft=8000;%采样频率 Fp1=1200; Fp2=3000; Fs1=1000; Fs2=3200; wp=[2*pi*Fp1/Ft,2*pi*Fp2/Ft]; ws=[2*pi*Fs1/Ft,2*pi*Fs2/Ft];%求出待设计的模拟滤波器的边界频率 [N,wn]=ellipord(wp,ws,Rp,Rs,'s'); [b,a]=ellip(N,Rp,Rs,wn,'s'); [bz,az]=bilinear(b,a,0.5); figure(7);%绘图 [h,w]=freqz(bz,az);%低通滤波器特性 title('IIR高通滤波器'); plot(w*fs/(2*pi),abs(h)); grid on %滤波 z=filter(bz,az,s); %滤波后画出其频谱，波形 Z=fft(z); %滤波后的信号频谱 figure(8) subplot(211); plot(z); title('带通滤波后的信号波形'); xlabel('时间'); ylabel('幅度'); grid on; subplot(212); plot(abs(fftshift(Z))); title('带通滤波后信号的频谱'); xlabel('频率'); ylabel('幅度'); axis([0 700000 0 20]); grid on; sound(z,fs); %回放滤波后的信号 ================运行后图像如下================