Search code examples
matlabsignal-processingcurve-fittingleast-squaresequalizer

Linear regression -- Stuck in model comparison in Matlab after estimation?


I want to determine how well the estimated model fits to the future new data. To do this, prediction error plot is often used. Basically, I want to compare the measured output and the model output. I am using the Least Mean Square algorithm as the equalization technique. Can somebody please help what is the proper way to plot the comparison between the model and the measured data? If the estimates are close to true, then the curves should be very close to each other. Below is the code. u is the input to the equalizer, x is the noisy received signal, y is the output of the equalizer, w is the equalizer weights. Should the graph be plotted using x and y*w? But x is noisy. I am confused since the measured output x is noisy and the model output y*w is noise-free.

%% Channel and noise level
h = [0.9 0.3 -0.1]; % Channel
SNRr = 10;              % Noise Level

%% Input/Output data
N = 1000;               % Number of samples
Bits = 2;               % Number of bits for modulation (2-bit for Binary modulation)
data = randi([0 1],1,N);        % Random signal
d = real(pskmod(data,Bits));    % BPSK Modulated signal (desired/output)
r = filter(h,1,d);              % Signal after passing through channel
x = awgn(r, SNRr);              % Noisy Signal after channel (given/input)

%% LMS parameters
epoch = 10;        % Number of epochs (training repetation)
eta = 1e-3;         % Learning rate / step size
order=10;           % Order of the equalizer

U = zeros(1,order); % Input frame
W = zeros(1,order); % Initial Weigths



%% Algorithm
for k = 1 : epoch
    for n = 1 : N
        U(1,2:end) = U(1,1:end-1);  % Sliding window
        U(1,1) = x(n);              % Present Input

        y = (W)*U';             % Calculating output of LMS
        e = d(n) - y;           % Instantaneous error 
        W = W +  eta * e * U ;  % Weight update rule of LMS
        J(k,n) = e * e';        % Instantaneous square error
    end
end

Solution

  • Lets start step by step:

    1. First of all when using some fitting method it is a good practice to use RMS error . To get this we have to find error between input and output. As I understood x is an input for our model and y is an output. Furthermore you already calculated error between them. But you used it in loop without saving. Lets modify your code:

      %% Algorithm
      for k = 1 : epoch
          for n = 1 : N
              U(1,2:end) = U(1,1:end-1);  % Sliding window
              U(1,1) = x(n);              % Present Input
      
              y(n) = (W)*U';             % Calculating output of LMS
              e(n) = x(n) - y(n);           % Instantaneous error 
              W = W +  eta * e(n) * U ;  % Weight update rule of LMS
              J(k,n) = e(n) * (e(n))';        % Instantaneous square error
          end
      end
      

      Now e consists of errors at the last epoch. So we can use something like this:

      rms(e)
      

      Also I'd like to compare results using mean error and standard deviation:

      mean(e)
      std(e)
      

      And some visualization:

      histogram(e)
      

      enter image description here

    2. Second moment: we can't use compare function just for vectors! You can use it for dynamic system models. For it you have to made some workaround about using this method as dynamic model. But we can use some functions as goodnessOfFit for example. If you want something like error at each step that consider all previous points of data then make some math workaround - calculate it at each point using [1:currentNumber].

    3. About using LMS method. There are built-in function calculating LMS. Lets try to use it for your data sets:

      alg = lms(0.001);
      eqobj = lineareq(10,alg);
      y1 = equalize(eqobj,x);
      

      And lets see at the result:

      plot(x)
      hold on
      plot(y1)
      

      enter image description here There are a lot of examples of such implementation of this function: look here for example.

    I hope this was helpful for you!