mydspexp1_bsdu.m

用matlab写的维纳滤波

clear all;
clc;
close all;
    h0 = [-2:0.01:4];
    h1 = [-4:0.01:2];
    [y,x] = meshgrid(h1,h0);
    J = 0.55 + x.^2 + y.^2 + 2* x .* y * cos(pi/8) ...
       - sqrt(2) * x * cos(pi/10) - sqrt(2) * y * cos(9*pi/40);
    Figure(1);
    mech(h0, h1, J);
    grid on;
                                                           %误差性能曲面
    Figure(2);
    contour(x, y, J,150);  hold on;                        %等值线
    plot(1.2 , h1);hold on;
    plot(h0 , -0.571);hold on;
    xlabel('h0');ylabel('h1');
%=========================================================================
%以上为等值线
%=========================================================================
clear all;
clc;
k = 0 : 1;
ItNum = 2000;                                               %递归次数
DelTa = 0.4;                                                %迭代步长
rxx = cos(2 * pi * k / 16);                                 %Rxx
ryx = cos(2 * pi * k / 16 + pi /10)/sqrt(2);                
Ryy0 = 0.55;                                                %Ryy(0)
Rxx = [rxx(1) rxx(2);rxx(2) rxx(1)];                        %Rxx自相关矩阵
Ryx = [ryx(1) ; ryx(2)];                                    %Ryx自相关矩阵
n = 0 : (ItNum - 1);
n0 = sin(2 * pi * n/16 + pi /10);                           %信号s
x = sqrt(2) * sin(2 * pi * n/16);                           %信号x
AVG = 100;                                                  %平均次数
%==========================================================================
%最陡下降法
%==========================================================================
% 
 Vg = zeros(2,1);                                            %梯度矩阵
 H = zeros(2,ItNum);
 H(:,1) = [3 ; -4];
 HalfDelTa = DelTa / 2.0;
 for i = 1 : (ItNum-1)
     Vg = Rxx * H(:,i) - Ryx;
     H(:,i+1) = H(:,i) - HalfDelTa * Vg;
 end
 plot(H(1,:),H(2,:),'r'); 
 hold on;
%==========================================================================    
% LMS算法
%==========================================================================
% Vg1 = zeros(2 , 1);                                       %临时梯度变量
error = zeros(1 , ItNum);                                   %存储噪声 
errorSQRTStore = zeros(1 , ItNum);                          %噪声平均
HStore = zeros(2 , ItNum);                                  %H滤波器系数矩阵
H1 = zeros(2 , ItNum);                                      %H计算矩阵
for i = 1 : AVG
    H1(:,1) = [3 ; -4];error(1) = 0;
    Noise = randn(1 , ItNum)/sqrt(20);                       %噪声
    y = n0 + Noise;                                          %y
    for j = 1 : (ItNum-1)
        error(j + 1)  = y(j + 1) - H1(: , j)' * [x(j + 1);x(j)];
        H1(:,j + 1) = H1( : , j) + DelTa * error(j + 1) * [x(j + 1) ; x(j)];
    end
    HStore  = HStore + H1;
    errorSQRTStore = errorSQRTStore + error .^2;
end
 plot(H1(1 , :),H1(2 , :));                                 %一次实现
 hold on;
 HStore = HStore/AVG;
 errorSQRTStore = errorSQRTStore/AVG;                       %平均
 plot(HStore(1 , :) , HStore(2 , :),'g');                   %100次平均
 figure(2);
 subplot(311);
 plot([1:ItNum] , error);                                   %Noise
 title('e(n)');
 subplot(312);
 plot([1:ItNum] , error .^ 2);                              %e(n)平方对比
 title('J(n)');
 subplot(313);
 plot([1 : ItNum] , errorSQRTStore);                        %e(n)平均
 title('J(n)平均');