fs=22050;
f=fs*(0:511)/1024;
while 1
K = menu('语音信号数字滤波器','语音采集','加噪声并回放','Butterworth滤波','回放语音','退出程序');
% guide
if (K==1)
x1=wavread('E:\Matlab滤波器\Windows XP 启动.wav'); %读取语音信号的数据，赋给变量x1
% x1=wavread('E:\Matlab滤波器\recycle.wav'); %读取语音信号的数据，赋给变量x1
% x1=wavread('E:\Matlab滤波器\1.wav'); %读取语音信号的数据，赋给变量x1
% x1=wavread('E:\Matlab滤波器\录音.wav'); %读取语音信号的数据，赋给变量x1
sound(x1,22050); %播放语音信号
F0=fft(x1,1024); %1024点fft变换
figure(1);
subplot(2,1,1);
plot(x1); %画出滤波后的时域图
title('原始信号的时域');
subplot(2,1,2);
plot(f,abs(F0(1:512))); %origion频谱图
title('原始信号的频谱')
xlabel('Hz');
ylabel('幅值');
elseif K==2 %加噪声
Y = AWGN(x1,20,0);
sound(Y,22050); %播放语音信号
figure(2);
subplot(2,1,1);
plot(Y);%加噪声高的时域
title('带噪声的时域');
y2=fft(Y,1024); %滤波前fft变换
subplot(2,1,2);
plot(f,abs(y2(1:512))); %画出滤波前的频谱图
title('带噪声的频谱')
xlabel('Hz');
ylabel('幅值');
elseif K==3
% 构造butterworth滤波器
wp=0.15*pi;
ws=0.1*pi;
Rp=1;
Rs=15;
Fs=22050;
Ts=1/Fs;
wp1=2/Ts*tan(wp/2); %将模拟指标转换成数字指标
ws1=2/Ts*tan(ws/2);
[N,Wn]=buttord(wp1,ws1,Rp,Rs,'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(3)
plot(W*Fs/(2*pi),abs(H))
grid
xlabel('频率／Hz')
ylabel('频率响应幅度')
title('Butterworth')
% 开始滤波
f1=filter(bz,az,Y);%滤波
F1=fft(f1,1024); %滤波后fft变换
figure(4)
subplot(2,1,1);
plot(f1); %画出滤波后的时域图
title('滤波后的时域');
subplot(2,1,2);
plot(f,abs(F1(1:512))); %画出滤波后的频谱图
title('滤波后的频谱')
xlabel('Hz');
ylabel('幅值');
elseif (K==4)
sound(f1,22050); %回放语音信号
elseif (K==5)
break;
end
end
% 参考资料
%一、 BUTTORD Butterworth filter order selection.
% [N, Wn] = BUTTORD(Wp, Ws, Rp, Rs) returns the order N of the lowest
% order digital Butterworth filter that loses no more than Rp dB in
% the passband and has at least Rs dB of attenuation in the stopband.
% Wp and Ws are the passband and stopband edge frequencies, normalized
% from 0 to 1 (where 1 corresponds to pi radians/sample). For example,
% Lowpass: Wp = .1, Ws = .2
% Highpass: Wp = .2, Ws = .1
% Bandpass: Wp = [.2 .7], Ws = [.1 .8]
% Bandstop: Wp = [.1 .8], Ws = [.2 .7]
% BUTTORD also returns Wn, the Butterworth natural frequency (or,
% the "3 dB frequency") to use with BUTTER to achieve the specifications.
%
% [N, Wn] = BUTTORD(Wp, Ws, Rp, Rs, 's') does the computation for an
% analog filter, in which case Wp and Ws are in radians/second.
%
% When Rp is chosen as 3 dB, the Wn in BUTTER is equal to Wp in BUTTORD.
% 二、AWGN Add white Gaussian noise to a signal.
% Y = AWGN(X,SNR) adds white Gaussian noise to X. The SNR is in dB.
% The power of X is assumed to be 0 dBW. If X is complex, then
% AWGN adds complex noise.
%
% Y = AWGN(X,SNR,SIGPOWER) when SIGPOWER is numeric, it represents
% the signal power in dBW. When SIGPOWER is 'measured', AWGN measures
% the signal power before adding noise.
%
% Y = AWGN(X,SNR,SIGPOWER,STATE) resets the state of RANDN to STATE.
%
% Y = AWGN(..., POWERTYPE) specifies the units of SNR and SIGPOWER.
% POWERTYPE can be 'db' or 'linear'. If POWERTYPE is 'db', then SNR
% is measured in dB and SIGPOWER is measured in dBW. If POWERTYPE is
% 'linear', then SNR is measured as a ratio and SIGPOWER is measured
% in Watts.
% 三、BUTTAP Butterworth analog lowpass filter prototype.
% [Z,P,K] = BUTTAP(N) returns the zeros, poles, and gain
% for an N-th order normalized prototype Butterworth analog
% lowpass filter. The resulting filter has N poles around
% the unit circle in the left half plane, and no zeros.
% 四、ZP2TF Zero-pole to transfer function conversion.
% [NUM,DEN] = ZP2TF(Z,P,K) forms the transfer function:
%
% NUM(s)
% H(s) = --------
% DEN(s)
%
% given a set of zero locations in vector Z, a set of pole locations
% in vector P, and a gain in scalar K. Vectors NUM and DEN are
% returned with numerator and denominator coefficients in descending
% powers of s.
% 五、 LP2LP Lowpass to lowpass analog filter transformation.
% [NUMT,DENT] = LP2LP(NUM,DEN,Wo) transforms the lowpass filter
% prototype NUM(s)/DEN(s) with unity cutoff frequency of 1 rad/sec
% to a lowpass filter with cutoff frequency Wo (rad/sec).
% [AT,BT,CT,DT] = LP2LP(A,B,C,D,Wo) does the same when the
% filter is described in state-space form.
% 六、FFT Discrete Fourier transform.
% FFT(X) is the discrete Fourier transform (DFT) of vector X. For
% matrices, the FFT operation is applied to each column. For N-D
% arrays, the FFT operation operates on the first non-singleton
% dimension.
%
% FFT(X,N) is the N-point FFT, padded with zeros if X has less
% than N points and truncated if it has more.
%
% FFT(X,[],DIM) or FFT(X,N,DIM) applies the FFT operation across the
% dimension DIM.
%
% For length N input vector x, the DFT is a length N vector X,
% with elements
% N
% X(k) = sum x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.
% n=1
% The inverse DFT (computed by IFFT) is given by
% N
% x(n) = (1/N) sum X(k)*exp( j*2*pi*(k-1)*(n-1)/N), 1 <= n <= N.
% k=1
评论0