lms.m

应用LMS算法去除噪声

clc;clear all;close all;
%----------纯净信号是8ASK信号，载频100KHz，波特率1KHz，采样率600K
x = dmod([0,3,1,3,0,2],100e3,1e3,600e3,'ask',8)'; 
%----------加入指定信噪比的噪声-------------------------------------
snr = 3;
N = length(x);
ps = sum(x.^2)/N;
sigma_noise = abs(sqrt(ps/(10^(snr/10))));
nn = randn(1,N);
enn = sum(nn)/N;                          
nn = nn-enn;                              
noise = sigma_noise*nn';
s = x+noise;    % s是含噪声的信号
%--------------------去噪之前的信噪比及均方差------------------------
snrs = 20*log10(norm(x)/norm(s-x));
mmses = mmse(s-x);
%--------------------参数的初始化------------------------------------
u = 0.002;         % 步长
M = 32;            % 滤波器阶数
d = x;
itr = length(x);
en = zeros(itr,1); % 误差
wn = zeros(M,itr); % 每行代表一个加权参量,每列代表一次迭代
%------------------LMS算法------------------------------------------
tic;
for k = M:itr      % 第k次迭代
    x_tap = x(k:-1:k-M+1);
    en(k) = d(k)-wn(:,k-1)'*x_tap;
    wn(:,k) = wn(:,k-1)+2*u*en(k,1)*x_tap;
end
y1 = conv(wn(:,length(x)),s);
y = max(abs(x))/max(abs(y1))*y1(1:length(s)); % y是滤波后的信号
toc;
%-------------------去噪之后的信噪比及均方差----------------------------
snry = 20*log10(norm(x)/norm(y-x));
mmsey = mmse(y-x);
snradd=snry-snrs    % 信噪比增量   
mmseadd=mmses/mmsey % 均方差增加的倍数 
%--------------------波形分析-------------------------------------------
figure,
subplot(3,1,1),plot(x);title('纯净信号');
subplot(3,1,2),plot(s);title('含噪声信号');
subplot(3,1,3),plot(y);title('去噪后信号');
figure,
plot(en);xlabel('采样点');ylabel('误差值');title('误差en');