THMMY.gr

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
  
  % add new value in buffer
  u = shift( u, 1 );
  u(1) = u_n;
  
  
  
  % your code here 
  % refresh vector u
  
  % apply w to u
  y_n = w'*u;
  e = d_n - y_n;
  % your code here 
  % update w for next iteration
  w = w + mu * e * u;
  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
  
  u(end+1,1) = u_n;
  
  % apply w to u
  y_n = w'*u(end:-1:end-m+1);
  
  e = d_n - y_n;
  % update w for next iteration
  w = w + mu* e * u(end:-1:end-m+1);
  
  end

% Step1: split output to top and bottom
function [yt, yb] = step1( x, w )

    % length of x signal
    n = length(x);

    % Check if n is a power of 2
    if mod(n,2) ~= 0
        disp('The signal x must have a length equal to a power of 2. \n')
        return
    else
        % k ranges from 0 to n - 1  
        k = (0:n-1)';

        % Split output top(yt) and bottom(yb) initiating both parts with 0s
        yb = zeros(n/2,1);
        yt = zeros(n/2,1);

        % Construct the 2-part output
        for j=0:n/2-1
            yt(j+1) = sum(w(n,j*k) .* x(k+1));      
        end 
    
        for j=0:n/2-1
            yb(j+1) = sum(w(n,(j+n/2)*k) .* x(k+1));
        end
        
        y = [yt ; yb];
    end
    
end

function [ybe, ybo, yte, yto] = step2( x, w )
    
    % split each output to even and odd
    n = length(x);

    % k ranges from 0 to n/2 -1  
    k = (0:n/2-1)';

    % Initiation of 4-part output with 0s 
    ybe = zeros(n/2,1);
    ybo = zeros(n/2,1);
    yte = zeros(n/2,1);
    yto = zeros(n/2,1);

    for j=0:n/2-1
        yte(j+1) = sum(w(n/2,j*k) .* x(2*k+1));
        yto(j+1) = w(n,j)*sum(w(n/2,j*k) .* x(2*k+2));
    end
    
    for j=0:n/2-1
        ybe(j+1) = sum(w(n/2,(j+n/2)*k) .* x(2*k+1));
        ybo(j+1) = w(n,j+n/2)*sum(w(n/2,(j+n/2)*k) .* x(2*k+2));
    end

    y = [yte + yto ; ybe + ybo];

end


function [xk1, xk2] = step3( x, w )
    
    % apply w identity and finish the proof with
    % y = [ xk1+xk2; xk1-xk2 ];
    n = length(x);
  
    % Initiation of output
    xk1 = zeros(n/2,1);
    xk2 = zeros(n/2,1);
  
    % k ranges from 0 to n/2 -1  
    k = (0:n/2-1)';

    for j=0:n/2-1
        xk1(j+1) = sum(w(n/2,j*k) .* x(2*k+1));
        xk2(j+1) = w(n,j)*sum(w(n/2,j*k) .* x(2*k+2));
    end

    y = [xk1+xk2; xk1-xk2];

end

function [y, flop] = rec_fast_fourier(x)

    % Initiate number of flops
    flop = 0;

    % n = length of signal to be transformed 
    n = length(x);

    % root of unity function handle
    w = @(n,k) exp(-2*k*pi*1i/n);

    % if n = 1 then the fast fourier transform cannot be broken into smaller parts
    % so its output is going to be an element of the original signal x.
    if n==1
        y = x;
    else    
        k = (0:n/2-1)';
        [yt,flopt] = rec_fast_fourier(x(2*k+1));
        [yb,flopb] = rec_fast_fourier(x(2*k+2));
        yb = w(n,k).*yb;
        y = [yt+yb ; yt-yb];
    
        flop = flopt + flopb + 6*n/2 + 2*n ;
    end
    
end