THMMY.gr

function [R, p] = getCorr(x, y, m)
 
  n = length(x);
 
  x = x(:);
  y = y(:);
   
  U = toeplitz( [x; zeros(m-1,1)], [x(1) zeros(1,m-1)] );
 
  y = [y; zeros(m-1,1)];
 
 
  R = U'*U / n;
  p = U'*y / n;
 
end

function [ybe, ybo, yte, yto] = step2( x, w )
% split each output to even and odd
 
  n = length(x);
 
  k = (0:n/2-1)';
 
  ybe = zeros(n/2,1);
  ybo = zeros(n/2,1);
  yte = zeros(n/2,1);
  yto = zeros(n/2,1);
 
  for m = 0:(n/2-1)
      yt(m+1) = x(m+1)*exp(-(2*pi*1i*m*k(m+1))/n);
  end
 
  for m = (n/2):(n-1)
      yb(m+1-n/2) = x(m+1)*exp(-(2*pi*1i*m*k(m+1-n/2))/n);
  end
 
  ybo = yb(:, 1:2:end);
  ybe = yb(:, 2:2:end);
  yto = yt(:, 1:2:end);
  yte = yt(:, 2:2:end);
 
end

function [yt, yb] = test( x )
% split output to top and bottom
 x=rand(16,1);
 %% all the initializations
  n = length(x);
  yt = zeros(n/2,1);
  yb = zeros(n/2,1);
  k = (0:n-1)';

  yt_after = zeros(n/2,1);
  yb_after = zeros(n/2,1);
  
 %% split in bottom and top
  yt = x((n/2):n); % y top
  yb = x(1:(n/2)); % y bottom
 
 %% do the dft for each part
  for ( k =1:1:n/2)
    s = 0;
    
    for (index=1:1:n/2)

        s = yb(index)*w(index,k) + s ;
    end
    
    yb_after(k) = s;
  end


  for ( k =1:1:n/2)
    s = 0;
    for ( index=1:1:n/2)
      s = yt(index)*w(index,k) + s ;
    end  
    yt_after(k) = s;
  end

  % save the values as outputs
    yb = yb_after;
    yt = yt_after;
  
    y =[ yb;yt];

    e=norm( y - fft(x));

   end

function [w]=w(n,k)

% this is my function handle w
% this is about the exp
% inputs : length of the time_samples n
% : length of the frequency_samples k

% what to do : only the exponential part
a = -i*2*pi*k*n;
b = 16; % should be 16 also the length of the given x
w = exp(a/16);
end

function [R] = solve( a, b, c, sigma_v1_sq, sigma_v2_sq )
  R = zeros(100,100);
  N=100;
  %Construct v1
  v1 = sqrt(sigma_v1_sq) * ones(N, 1);
  u = v1;
  %% Generating input signal
    u(2) =u(2) -a(1)*u(1);
    u(3) =u(3) -a(1)*u(2)-a(2)*u(1);
    for n = 4:N
        u(n)=u(n) -a(1)*u(n-1)-a(2)*u(n-2)-a(3)*u(n-3)+v1(n);
    end
  % Calculate autocorrelation (synarthsh apo dokimastiko quiz 2019)
  R = getCorr(u,N)
end  

function [R] = getCorr(x,m)  
n = length(x);
x = x(:);
U = toeplitz( [x; zeros(m-1,1)], [x(1) zeros(1,m-1)] );
R = U'*U/n;
end

function y_n = myLMS(u_n, d_n, m, mu)

  persistent w u
 
  if isempty(w)
    w = zeros(m,1);
    u = zeros(m,1);
  end
 
  % your code here
  % refresh vector u
 
  u = [u_n; u(1:m-1)]
 
  % apply w to u
  y_n = w'*u;
 
  % your code here
  % update w for next iteration
  e = d_n - y_n;
  w = w + mu * u * e;
 
  end
 

function [y,w_out] = rls(u, d, m, lambda, delta)
% RLS adaptation
%   
 
  %persistent w
  %persistent u_stream
  %persistent P
 
    u = u(:);
    d = d(:);
    n = length(u);

    P = delta * eye(m);
    w_out = zeros(m, 1);
    y = zeros(size(u));
    e = zeros(size(u));

    for i = (m + 1):n
        ubar = u(i:-1:i - m + 1);
        y(i) = w_out.' * ubar;
        k = 1 / lambda * P * ubar / (1 + 1 / lambda * ubar' * P * ubar);
        e(i) = d(i) - y(i);
        w_out = w_out + k * e(i);
        P = 1 / lambda * P - 1 / lambda * k * ubar.' * P;
    end
 
end

function [R, p, w_o] = wiener( a, b, ~, sigma_v1, ~ )

    %% Reformat args
    a = [1; a(:)];


    %% Υπολογισμός πίνακα αυτοσυσχέτισης
    A = [...
        a(1) a(2) a(3) a(4);...
        a(2) a(1)+a(3) a(4) 0;...
        a(3) a(2)+a(4) a(1) 0;...
        a(4) a(3) a(2) a(1);...
    ];
   
    B = [sigma_v1; 0; 0; 0];
   
    r = A \ B;
    R = toeplitz(r);    % πίνακας αυτοσυσχέτισης
   
    %% Υπολογισμός διανύσματος ετεροσυσχέτισης
    p = R * (-b');
   
    %% Υπολογισμός Wiener Coefficients
    w_o = R \ p;
 
end

function [y,w_out] = lms(u, d, m, mu)
% LMS - Least Mean Square adaptation
%  
 
  persistent w
  persistent u_stream
  
  % initialize w
  if isempty(w)
    w = zeros( m, 1 );
  end
  
  u_stream(end+1, 1) = u;
 
  if length( u_stream ) < m                
    y = 0;  
  else                                  
    % compute output
    y = w' * u_stream(end:-1:end-m+1);  
    % get error
    e = d - y;
    
    if length( u_stream ) > m    
      % update weights
      w = w + mu * e * u_stream( end:-1:end-m+1 );      
    end    
  end
  
  if nargout > 1
    w_out = w;
  end

end

[quote author=annaTZ link=topic=74039.msg1297525#msg1297525 date=1560251202]
Μπορει κανεις να μου εξηγησει τι κανω λαθος στην 4η ασκηση? Παιρνω 0.25, δηλαδη ο wiener βγαινει σωστος αλλα τα υπολοιπα οχι...
Δεν μπορω να καταλαβω γιατι, τα κανω οπως ο φλωρος στην απαντηση που ειχε δωσει σε προηγουμενη ασκηση

function [y,w_out] = rls(u, d, m, lambda, delta)
% LMS - RLS adaptation
%   
 
    persistent w u_stream P
    if isempty( u_stream )
       
        P = delta * eye( m );
        w = zeros( m, 1 );
        u_stream = [ u; zeros( (m - 1),1 ) ];
       
    else
       
        u_stream = [ u; u_stream( 1 : ( m - 1 ) ) ];
       
    end   
   
    %% LMS Part
    y = w' * u_stream;      % filter input
    e = d - y;            % compute error
   
    %% RLS Part
    k = 1 / lambda * P * u_stream / (1 + 1 / lambda * u_stream' * P * u_stream);
    P = 1 / lambda * P - ( 1 / lambda ) * k * u_stream.' * P;
    w = w + k * e;          % get next Wiener approximation
   
    w_out = w;
 
end

function [R,p,w_o] = wiener( a, b, c, sigma_v1, sigma_v2 )
 
    a = [1; a(:)];

    A = [a(1) a(2) a(3) a(4); a(2) a(1)+a(3) a(4) 0; a(3) a(2)+a(4) a(1) 0; a(4) a(3) a(2) a(1)];
   
    B = [sigma_v1; 0; 0; 0];
   
    r = A \ B;
    R = toeplitz(r);
   
    p = R * (-b');
   
    w_o = R \ p;
 
end


function [y,w_out] = lms(u, d, m, mu)
% LMS - Least Mean Square adaptation
%   
 
  persistent w
  persistent u_stream
 
 
  % initialize
  if isempty(w)
    w = zeros( m, 1 );
    u_stream = zeros(m,1);
  end
 
  u_stream(end+1, 1) = u;
                           
   
    % compute output
    y = w' * u_stream(end:-1:end-m+1);
   
    % get error
    e = d - y;
   
                     
   
      % update weights
      w = w + mu * e * u_stream( end:-1:end-m+1 );

 

  if nargout > 1
    w_out = w;
  end
 
 
end


function [y,w_out] = nlms(u, d, m, mu, alpha)
% LMS - Least Mean Square adaptation
%   
 
  persistent w
  persistent u_stream
 
 
  % initialize
  if isempty(w)
    w = zeros( m, 1 );
    u_stream = zeros(m,1);
  end
 
  u_stream(end+1, 1) = u;

 
                       
   
    % compute output
    y = w' * u_stream(end:-1:end-m+1);
   
    % get error
    e = d - y;

      % update weights
      w = w + mu * e * u_stream( end:-1:end-m+1 )/ (alpha + u_stream( end:-1:end-m+1 )' * u_stream( end:-1:end-m+1 ));
     

   


  if nargout > 1
    w_out = w;
  end
 
 
end



function [y,w_out] = rls(u, d, m, lambda, delta)
% LMS - RLS adaptation
%   
   
  persistent w
  persistent u_stream
  persistent P
 
 
  % initialize
  if isempty(w)
    w = zeros( m, 1 );
    u_stream = zeros(m,1);
  P = (1 / delta) * eye(m, m);
     
  end
 
  u_stream(end+1, 1) = u;

                         
   
    % compute output
    y = w' * u_stream(end:-1:end-m+1);
   
    % get error
    e = d - y;
   
 
      %compute k
      k = ((lambda^-1) * P * u_stream( end:-1:end-m+1 ) / (1 + (lambda^-1) * u_stream( end:-1:end-m+1 )'* P * u_stream( end:-1:end-m+1 )));

                               
                                 
      % update weights
      w = w + k * e;
     
      %compute P
      P = (lambda^-1) * P - (lambda^-1) * k * u_stream( end:-1:end-m+1 )' * P;
     


  if nargout > 1
    w_out = w;
  end
 
 
end