Search code examples
matlabsignal-processingconvolutionestimationdeconvolution

Matlab : Convolution and deconvolution results weird


Data x is input to an autoregreesive model (AR) model. The output of the AR model is corrupted with Additive White Gaussian Noise at SNR = 30 dB. The observations are denoted by noisy_y.

Let there be close estimates h_hat of the AR model (these are obtained from Least Squares estimation). I want to see how close the input obtained from deconvolution with h_hat and the measurements is to the known x.

  • My confusion is which variable to use for deconvolution -- clean y or noisy y?

Upon deconvolution, I should get x_hat. I am not sure if the correct way to perform deconvolution is using the noisy_y or using the y before adding noise. I have used the following code.

  • Can somebody please help in what is the correct method to plot x and x_hat.

Below is the plot of x vs x_hat. As can be seen, that these do not match. Where is my understand wrong? Please help.

image

The code is:

clear all
    N = 200; %number of data points
    a1=0.1650;
    b1=-0.850;
    h = [1 a1 b1]; %true coefficients

    x = rand(1,N);
    %%AR model
    y = filter(1,h,x); %transmitted signal through AR channel
    noisy_y = awgn(y,30,'measured');
    hat_h= [1 0.133 0.653];
    x_hat = filter(hat_h,1,noisy_y); %deconvolution
    plot(1:50,x(1:50),'b');
    hold on;
    plot(1:50,x_hat(1:50),'-.rd');

Solution

  • A first issue is that the coefficients h of your AR model correspond to an unstable system since one of its poles is located outside the unit circle:

    >> abs(roots(h))
    ans =
    
       1.00814
       0.84314
    

    Parameter estimation techniques are then quite likely to fail to converge given a diverging input sequence. Indeed, looking at the stated hat_h = [1 0.133 0.653] it is pretty clear that the parameter estimation did not converge anywhere near the actual coefficients. In your specific case you did not provide the code illustrating how you obtained hat_h (other than specifying that it was "obtained from Least Squares estimation"), so it isn't possible to further comment on what went wrong with your estimation.

    That said, the standard formulation of Least Mean Squares (LMS) filters is given for an MA model. A common method for AR parameter estimation is to solve the Yule-Walker equations:

    hat_h = aryule(noisy_y - mean(noisy_y), length(h)-1);
    

    If we were to use this estimation method with the stable system defined by:

    h = [1 -a1 -b1];
    
    x = rand(1,N);
    %%AR model
    y = filter(1,h,x); %transmitted signal through AR channel
    noisy_y = awgn(y,30,'measured');
    hat_h = aryule(noisy_y - mean(noisy_y), length(h)-1);
    x_hat = filter(hat_h,1,noisy_y); %deconvolution
    

    The plot of x and x_hat would look like:

    enter image description here