MATLAB对语音信号加随机噪声及去噪程序

[y,fs,bits]=wavread(‘C:\Documents and Settings\Administrator\桌面\cuocuo.wav‘);

%  sound(y,fs)      % 回放语音信号

n=length(y)  %选取变换的点数 

y_p=fft(y,n);      %对n点进行傅里叶变换到频域

f=fs*(0:n/2-1)/n;   % 对应点的频率

figure(1)

subplot(2,1,1);

plot(y);                    %语音信号的时域波形图

title(‘原始语音信号采样后时域波形‘);

xlabel(‘时间轴‘)

ylabel(‘幅值 A‘)

subplot(2,1,2);

plot(f,abs(y_p(1:n/2)));     %语音信号的频谱图

title(‘原始语音信号采样后频谱图‘);

xlabel(‘频率Hz‘);

ylabel(‘频率幅值‘);

%对音频信号产生噪声

  L=length(y)        %计算音频信号的长度

  noise=0.1*randn(L,2);  %产生等长度的随机噪声信号(这里的噪声的大小取决于随机函数的幅度倍数）

  y_z=y+noise;        %将两个信号叠加成一个新的信号——加噪声处理   

  %sound(y_z,fs)

%对加噪后的语音信号进行分析

n=length(y);  %选取变换的点数 

y_zp=fft(y_z,n);      %对n点进行傅里叶变换到频域

f=fs*(0:n/2-1)/n;   % 对应点的频率

figure(2)

subplot(2,1,1);

plot(y_z);                    %加噪语音信号的时域波形图

title(‘加噪语音信号时域波形‘);

xlabel(‘时间轴‘)

ylabel(‘幅值 A‘)

subplot(2,1,2);

plot(f,abs(y_zp(1:n/2)));     %加噪语音信号的频谱图

title(‘加噪语音信号频谱图‘);

xlabel(‘频率Hz‘);

ylabel(‘频率幅值‘);

对加噪的语音信号进行去噪程序如下：

fp=1500;fc=1700;As=100;Ap=1;

(以上为低通滤波器的性能指标）

wc=2*pi*fc/fs; wp=2*pi*fp/fs;

wdel=wc-wp;

beta=0.112*(As-8.7);

N=ceil((As-8)/2.285/wdel);

wn= kaiser(N+1,beta);

ws=(wp+wc)/2/pi;

b=fir1(N,ws,wn);

figure(3);

freqz(b,1);         

（此前为低通滤波器设计阶段）——接下来为去除噪声信号的程序——

x=fftfilt(b,y_z);

X=fft(x,n);

figure(4);

subplot(2,2,1);plot(f,abs(y_zp));

title(‘滤波前信号的频谱‘);

subplot(2,2,2);plot(f,abs(X));

title(‘滤波后信号频谱‘);

subplot(2,2,3);plot(y_z);

title(‘滤波前信号的波形‘)

subplot(2,2,4);plot(x);

title(‘滤波后信号的波形‘)

%sound(x,fs,bits)  %回放滤波后的音频

      设计滤波器：

器常用的方法有：脉冲响应不变法和双线性变换法。

数字信号中我们了解到，脉冲响应不变法好处是模拟角频率Ω和数字频率w呈现线性关系w=ΩT，但主要的缺点是它会产生频谱混叠现象，使数字滤波器的频率响应偏离模拟滤波器的频响特性。为了克服这一缺点我们常采用双线性变换法。

下面我们给出IIR数字滤波器和FIR数字滤波器的设计步骤：

一、IIR滤波器的设计步骤如下：

  （1）、确定数字低通滤波器的技术指标：通带边界频率、通带最大衰减、阻带截止频率、阻带最小衰减。

  （2）、将数字低通滤波器的技术指标转换成响应的模拟低通滤波器的技术指标。

  （3）、按照模拟低通滤波器的技术指标设计模拟低通滤波器。

  （4）、用双线性变换法将模拟滤波器系统函数转换成数字低通滤波器系统函数。

二、FIR滤波器的设计步骤如下：（重点介绍窗函数法）

  （1）、根据对阻带衰减以及过渡带的指标要求，选择窗函数类型（有：矩形窗、三角窗、汉宁窗、哈明窗、凯塞窗等），并估计窗口长度N。先按照阻带衰减选择窗函数类型。在保证阻带衰减满足要求的情况下，尽量选择主瓣集中的窗函数。

  （2）、构造希望逼近的频率响应函数。

  （3）、计算h(n)，有限长序列

  （4）、加窗得到设计结果。