wiener_1.m

来自「一维WIENER的方法.包括了滤波器的频谱图和幅度响应图.可直接使用」· M 代码 · 共 35 行

clear;

%信号的生成
maxlag=100;
w0=3000;
fs=1000;
t=0:1/fs:(1-1/fs);
r=rand(1,fs);
x=sin(w0*t+r);%源信号，要具平稳性，产生平稳信号见通信原理书 P18
y=x+0.8*randn(1,fs);%加噪后的信号

%wiener filter的生成
[rx,lags]=xcorr(y,maxlag,'biased');%观测信号的自相关函数
rx1=toeplitz(rx(101:end));%对称化自相关函数矩阵使之成为方阵，filter的阶数为101阶
rx2=xcorr(x,y,maxlag,'biased');%观测信号与期望信号的互相关函数
rx2=rx2(101:end);
h=inv(rx1)*rx2';%Wiener-Horf 方程
y0=filter(h,1,y);%加噪信号通过滤波器后的输出

%波形输出
s1=fft(x,1000);
s1=fftshift(s1);
s2=fft(y,1000);
s2=fftshift(s2);
s3=fft(y0,1000);
s3=fftshift(s3);
f=linspace(-fs/2,fs/2,1000);
figure(1)
subplot(1,3,1);stem(f,abs(s1));title('原信号频谱');axis tight%源信号的频谱
subplot(1,3,2);stem(f,abs(s2));title('加噪信号频谱');axis tight%加噪信号的频谱
subplot(1,3,3);stem(f,abs(s3));title('滤波器输出信号频谱');axis tight%滤波后信号的频谱

[h101 w101]=freqz(h,1,101*101);
figure(2)
stem(w101,abs(h101));title('winer滤波器的幅度响应')