2

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');
SKM
  • 959
  • 2
  • 19
  • 45

1 Answers1

1

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

SleuthEye
  • 14,379
  • 2
  • 32
  • 61
  • Thank you for your answer. I have 4 doubts, could you please clarify? (1)I was not aware of the unstable problem. If `h = [1, 0.195,-0.95];` then taking the `abs(root(h)) = 1.0770, 0.8820` Is this unstable as well? If the first value, 1.0770 > 1, then do we say that the root is outside the unit circle?(2) Error calculation: using the `h_hat` and `x_hat` I will simulate the AR model to generate `y_hat`. Since the estimates are obtained using `noisy_y`, should the mean square error be calculated between `y` and `y_hat` or would it be `noisy_y` and `y_hat` ? – SKM Oct 26 '17 at 19:27
  • (3) Should I sub-divide the data sequence into training and testing to estimate the system and calculate the `y_hat`? (4) When the input and the system parameters are unknown, I don't have desired signal / training sequence in the form of the pilot symbols of the input. In this case, would estimation and `y_hat` calculation be done using the entire sequence (like shown in this Question)? – SKM Oct 26 '17 at 19:27
  • 1
    (1) `h=[1,0.195,-0.95]` is unstable as well. It is particularly obvious if you plot a long `x` sequence. If any values in `abs(roots(h))` is greated than 1 then the root is outside the unit circle. (2) probably better as a different question where you expose exactly what you are trying to achieve with the comparison, but you should in all likelihood be comparing like things (e.g. `y` with `y_hat` and `noisy_y` with `noisy_y_hat`) – SleuthEye Oct 27 '17 at 00:45
  • 1
    (3) The first few elements of `y_hat` might not be so accurate until the system parameters have reasonable estimates; hence an initial training phase usually makes sense. (4) Blind estimation may take more samples to converge than a pilot aided one, but the exact length of the sequence would depend on the exact estimation method used and desired performance. – SleuthEye Oct 27 '17 at 00:45
  • Thank you for answering to all my comments. One last thing which I wanted to ask from you is that the unit circle condition is also applicable to moving average model? – SKM Oct 27 '17 at 02:21
  • 1
    No, stability conditions only requires that poles of the transfer function be located inside the unit circle. Moving average model doesn't have poles, but only zeros (hence they are always stable). – SleuthEye Oct 27 '17 at 02:44
  • Thank you for your answers in the comment. Sorry to bother you again. But I closely looked at one of your comments and just wanted to confirm if I understood correctly or not. In your comment you said tht like things should be compared such as `y_hat` and `y`. I will get `y_hat` using convolution between `h_hat` and `x_hat` which are the estimates obtained from `noisy_y` for which you have put up the correct code. Is this `y_hat`, the `noisy_yhat` because it is obtained from the estimates from the noisy observations? – SKM Oct 27 '17 at 06:39
  • I have posted a new question for this here https://stackoverflow.com/questions/46969316/matlab-comparison-after-deconvolution Shall be immensely grateful for your help if you can answer this as well. – SKM Oct 27 '17 at 07:05
  • I said "in all likelyhood" because it didn't know at the time how you would be calculating `y_hat` and assumed it might be `y_hat = filter(1,hat_h,x); noisy_y_hat = awgn(y_hat,30,'measured')`, but I see from [your follow-up question](https://stackoverflow.com/questions/46969316/matlab-comparison-after-deconvolution) that it isn't the case. – SleuthEye Oct 27 '17 at 12:15
  • So do you think I should compare `y_hat = filter(1,hat_h,hat_x)` with `y`? Maybe if you can answer the new question, it would be helpful since you already have a better picture of the problem.:) thank you for your time and effort. – SKM Oct 27 '17 at 17:48