lmsyz.m

自适应滤波算法的LMS算法的仿真验证程序是在FM信号上提取完成的

clear;
% close all;
ts=cputime;
clc;
sysorder=10   %滤波器阶数；
fc=2e4;       
fd=5000          %搬移频率；
[fn,fp]=uigetfile('*.wav','pick a wav-file');
[y,fs,bit]=wavread([fp,fn]);
for a=1:fs
yy(a,1)=y(a,1);
end
n=1:fs;
opt=(fc/fs)*2*pi/(max(max(yy)));
%  x=modulate(yy,fc,fs,'fm',opt);
x=cos(2*pi*fc*n'/fs+opt*cumsum(yy));  %信号调制；
carr=cos(2*pi*fd*n'/fs);
adata=x.*carr;            %信号搬移；
sum=zeros(1,fs);
      h=[0.0976 
          0.2873
          0.3360 
          0.2210 
          0.0964];        %权矢量期望值；
     x1=conv(h,x);
     d=x;      
     s=adata+x1(1:fs);  %主回路信号；
     w = zeros ( sysorder ,1 ) ;       %权矢量初始化；
   u=1/(s'*s/fs*sysorder)*0.01;        %步长；
   % 开始算法；
      for n = sysorder : fs
 U = d(n:-1:n-sysorder+1) ; 
 y(n)= w' * U; 
e(n) = s(n) - y(n) ; 
w = w + 2*u * U * e(n) ; 
      end
 %算法验证；
 for n = sysorder :fs 
U = d(n:-1:n-sysorder+1) ; 
y(n) = w'* U ; 
e(n) = s(n) - y(n) ; 
 end
 ee=e';
xcorr1=xcorr(x,s);
abs1=abs(xcorr1);
xcorr2=xcorr(x,ee);
abs2=abs(xcorr2);
xcorr3=xcorr(x,adata);
abs3=abs(xcorr3);
 figure(1)
subplot(311)
[px1,f1]=pwelch(yy,[],0,512,fs);
plot(f1,px1);
title('未调制信号的功率谱');
xlabel('频率(Hz)');
ylabel('功率谱');
subplot(312)
[px2,f2]=pwelch(x,[],0,512,fs);
plot(f2,px2);
title('信号调制后');
xlabel('频率(Hz)');
ylabel('功率谱');
subplot(313)
[px3,f3]=pwelch(adata,[],0,512,fs);
plot(f3,px3);
title('信号搬移后')
xlabel('频率(Hz)');
ylabel('功率');
figure(2)
subplot(311)
[px1,f1]=pwelch(s,[],0,512,fs);
plot(f1,px1);
title('信号叠加');
xlabel('频率(Hz)');
ylabel('功率谱');
subplot(312)
[px2,f2]=pwelch(d,[],0,512,fs);
plot(f2,px2);
title('调制后信号');
xlabel('频率(Hz)');
ylabel('功率谱(db)');
subplot(313)
[px3,f3]=pwelch(ee,[],0,512,fs);
plot(f3,px3);
title('误差');
xlabel('频率(Hz)');
ylabel('功率谱');
figure(3)
subplot(211);
plot(abs1);
title('叠加信号与调制后信号的相关性');
axis([0 5e4 0 6000]);
subplot(212);
plot(abs2);
title('调制信号与误差信号的相关特性');
axis([0 5e4 0 6000]);
figure(4)
plot(abs3);
title('调制信号与期望信号的相关特性');
axis([0 5e4 0 6000]);
figure(5)
hist(x);
figure(6)
hist(adata);
figure(7)
plot(h, 'k+');
hold on
plot(w, 'r*');
legend('Actual weights','Estimated weights');
title('Comparison of the actual weights and the estimated weights') ;
xlabel('权矢量阶数');
ylabel('权值');