MATLAB 数字信号处理 滤波器设计
帅小龙 2019-11-29 09:26:31 设计目的和内容
通过数字信号处理的大作业设计,巩固和运用数字信号处理课程中的理论知识和实验技能,掌握最基本的数字信号处理的理论和方法。本次设计的题目是:语音信号的处理与滤波。
具体来说就是对一段语音信号进行采样;画出采样后语音信号的时域波形和频谱图;然后给采样信号加以噪声,采用双线性变换设计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); %回放滤波后的信号
================运行后图像如下================