音频信号处理频谱混叠滤波

信号与系统综合训练项目报告1

班级电气14-5

姓名___王怡然史晓玉

学号14050405201405040514

指导教师李楠信号与系统综合训练项目一

声音信号的频谱分析与处理

1目的

1.加深对信号频域分析基本概念、基本理论和基本方法的理解和掌握

2.初步学会使用matlab以及音频工具进行声音信号的分析。

3.通过观察声音信号的频谱，加深对实际信号分析问题的认识和理解，学会运用fft等理论工具分析信号频谱。

4.使用低筒滤波器对音频信号进行处理，并听取处理后音频效果，与之前信号进行对比，得出结论。

5.利用matlab将该音频信号产生带回音的音频信号，消除带回音的音频信号，掌握系统以及信号处理技术在声音信号处理中的应用。

2原始信号的时域波形，FFT频谱分析

1.1原始声音打开

在MATLAB中，[x,fs,bits]=wavread('2014102419-15-33.wav');用于读取语音，sound(x,fs,bits);用于对声音的回放。

gequnn.wav

2.1相关程序

在MATLAB中，[y,fs,bits]=wavread('Blip',[N1N2]);用于读取语音，采样值放在向量y中，fs表示采样频率(Hz)，bits表示采样位数。[N1N2]表示读取从N1点到N2点的值（若只有一个N的点则表示读取前N点的采样值）。sound(x,fs,bits);用于对声音的回放。向量y则就代表了一个信号（也即一个复杂的“函数表达式”）也就是说可以像处理一个信号表达式一样处理这个声音信号。利用函数fft对信号行快速傅里叶变换,得到信号的频谱特性，从而加深对频谱特性的理解。n=length(y)，求出语音信号的长度。Y=fft(y,n)，快速傅里叶变换。grid添加网络。subplot把画图区域划分。

程序如下：

[x,fs,bits]=wavread('C:\\Users\\Administrator.USER01604261154\\Desktop\\g equnn.wav');

sound(x,fs,bits);

n=length(x);

figure(1);

subplot(3,1,1);

plot(x);grid on;

xlabel('时间');ylabel('幅值');title('原始信号时域波形');

subplot(3,1,2);

wx=fft(x);

plot(abs(wx));grid on;

xlabel('频率');ylabel('幅值');title('幅频特性');

subplot(3,1,3);

plot(angle(wx));grid on;

xlabel('频率');ylabel('相位');title('相频特性');

2.1.1频谱分析：

3白噪声的产生

白噪声就是频谱为常数，协方差函数在delay=0时不为0，在delay不为0时等于0，即样本点互不相关。当随机的从高斯分布中获取采样值时，采样点所组成的随机过程就是“高斯白噪声”。而MATLAB中函数y=randn(1,length(t))是一个产生length(t)个随机数的函数，即从高斯分布中获取样值并集合采样点的随机过程，所以它可以用来产生高斯白噪声。将高斯白噪声叠加到原始的音频信号上。

3.1加噪声程序：

[y,fs,nbits]=wavread('C:\\Users\\Administrator.USER01604261154\\Desktop\\ gequnn.wav');

n=length(y);

Noise=0.2*randn(n,2);

s=y+Noise;

sound(s,fs,nbits);

figure(1);

subplot(2,1,1);

plot(s);

title('加噪语音信号的波形')

grid;

S=fft(s);

subplot(2,1,2);

plot(abs(S));

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

grid;

3.1.1加噪声后的频谱分析

IIR滤波器原理

IIR是无限冲激响应滤波器，原理是进行卷积计算，即对输入信号进行卷积计算。冲激响应不变法是使数字滤波器在时域上模拟滤波器，但是它们的缺点是产生频率响应的混叠失真，这是由于从s平面到z平面是多值的映射关系所造成的。

双线性变换法是使数字滤波器的频率响应与模拟滤波器的频率响应相似的一种变换方法。为了克服多值映射这一缺点，我们首先把整个s平面压缩变换到某一中介的s1平面的一条横带里，再通过变换关系将此横带变换到整个z平面上去，这样就使得s平面与z平面是一一对应的关系，消除了多值变换性，也就消除了频谱混叠现象。

双线性法设计IIR数字滤波器的步骤：

1）将数字滤波器的频率指标{k}由Wk＝（2/T）*tan(wk),转换为模拟滤波器的频率指标{k}.2)由模拟滤波器的指标设计由H(s)转换为H(z)。

3.2去噪声程序

[y,fs,nbits]=wavread('C:\\Users\\Administrator.USER01604261154\\Desktop\\ gequnn.wav');

n=length(y);

Noise=randn(n,2);

s=y+0.2*Noise;

S=fft(s);

fs=22000;

Ts=1/fs;R1=10;

wp=2*pi*1000/fs;%通带截止频率

ws=2*pi*1200/fs;%阻带截止频率

Rp=1;%通带衰减

Rl=100;%阻带衰减

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

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

[N,Wn]=buttord(wp1,ws1,Rp,R1,'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('IIR低通滤波器')

f1=filter(bz,az,y);

figure(2);

subplot(2,1,1);

plot(y);%画出滤波前的时域图title('IIR低通滤波器滤波前的时域波形');

subplot(2,1,2)

plot(f1);%画出滤波后的时域图title('IIR低通滤波器滤波后的时域波形');

sound(f1,44100);%播放滤波后的信号F0=fft(f1,1024);

f=fs*(0:511)/1024;

y2=abs(fft(f1));

figure(4);

subplot(2,2,3);

plot(f1);

title('滤波后信号的时域图');

subplot(2,2,4);

plot(y2);

title('滤波后信号的频谱');

3.2.1去噪声频谱分析

3.3结果分析

将加噪声的频谱和滤波后的频谱做对比，发现滤波后的音乐信号发生了衰减，说明滤波器起到了滤波作用，同时通过幅频对比，可以看出滤波器滤掉了一部分频率范围内的信号。分别听原始声音信号和滤波后的声音信号，发现滤波后的声音信号单频噪声消失，同时原始信号强度稍有减弱。

4对原始信号加回音

4.1声音信号y[n]的产生

因为人耳能分辨的声音延迟至少0.1秒，因此，最小延迟不能小于0.1秒。在已有声音信号x的基础上产生带回声的声音信号，可以表达为在于原信号的基础上叠加其延时的分量。假设只有一个回声的情况下，可简化其模型为：y(n)=x(n)+ax(n-N),a为反射系数，N为延时时间。在此先延时3秒，即最小延迟量N=3*22000=66000。这里设N为66000，a=0.2

4.2利用矩阵置零产生x的延迟，以及得到y信号

为了保证图像的完整性，对读取的信号先延长76000个采样点，将原来信号延长66000个采样点，然后在后面补上10000个点。

4.3加回声程序：

[x,fs,nbits]=wavread('C:\\Users\\Administrator.USER01604261154\\Desktop\\ gequnn.wav');%把语音信号进行加载入Matlab仿真软件平台中。

N=length(x);%语音信号的长度。

x1=x(1:N);

x2=x(1:N);

x1=[x1,zeros(1,76000)];%zeros(1,76000)产生1行76000列全零矩阵加到x1后面。x2=[zeros(1,66000),0.2*x2,zeros(1,10000)];%N+66000+x=N+76000，得

x=10000

y=x1+x2;%加入回音的信号。

figure(2);

subplot(3,1,1);

plot(y(1:66000.+N));

title('含回声信号波形');

y1=fft(y);

subplot(3,1,2);

plot(abs(y1));

title('含回声信号幅值');

subplot(3,1,3);

plot(angle(y1));

title('含回声信号相位');

sound(y,fs);

4.4去回声程序：

[x,fs,nbits]=wavread('C:\\Users\\Administrator.USER01604261154\\Desktop\\ jiequ.wav');%把语音信号进行加载入Matlab仿真软件平台中。

N=length(x);%语音信号的长度。

x1=x(1:N);

x2=x(1:N);

x1=[x1,zeros(1,76000)];%zeros(1,76000)产生1行76000列全零矩阵加到x1后面。x2=[zeros(1,66000),0.2*x2,zeros(1,10000)];%N+66000+x=N+76000，得

x=10000

y=x1+x2;%加入回音的信号。

a=[1,zeros(1,9999),0.2];

b=1;

x1=filter(b,a,y);%滤波函数，b,a分别系统函数分子分母多项式的系数，结果放入y wavwrite(x1,'恢复后的信号')

figure(3);

subplot(3,1,1);

plot(x1);grid on;

xlabel('时间');ylabel('幅值');title('恢复后的信号时域波形');

subplot(3,1,2);

wx1=fft(x1);

plot(abs(wx1));grid on;

xlabel('频率');ylabel('幅值');title('幅频特性');

subplot(3,1,3);

plot(angle(wx1));grid on;

xlabel('频率');ylabel('相位');title('相频特性');

sound(x1,fs,nbits);

4.5结果分析

从时域上看，波形和滤波前大体相同，幅度上有微小变化。从频谱上看，幅

值略微减小，曲线总体趋势向下，相位部分地方跳动明显。