[Ch11. Parametric Signal Model (a.k.a. AR Model)] Better Frequency Resolution with smaller points when estimating Frequency of Sinusoidal

    When it comes to estimating the main frequency of sinusoidal like signals with finite length, the most intuitive method might be observing the largest component in the DFT/DTFT. Nonetheless, the resolution of DFT is limited by windowing effect. If a known structure of the signals is given, these method are not the best way to estimate the frequency.

    With the known structure of the signals, Professor Oppenheim introduced an alternative way to estimate frequency, Parametric Signal Modeling (a.k.a. Autoregressive model, AR Model)[1]. To apply this method, a structure of all-pole linear system, H(z), should be predefined. This all-pole system tries to represent the signal with relative small set of parameters (or coefficients). Then, we estimate the coefficients of the system with observed signals. Last, the DTFT of the system can be regarded as the power spectrum density (PSD) of the observed signals.

    There are two basic methods to estimate coefficients from a observed signal s[n] with  finite length. One is autocorrelation method and the other is covariance method. The main difference between these two methods is the definition of prediction error. The length of s[n] is M+ 1 and the order of H(z) is p.

  • Autocorrelation Method considers all M+p length of error (i.e. n ranges from zero to M+p).
    • Prediction error of the first p and the last p samples can be large due to zero samples.
    • Hamming window is applied to reduced the influence, but the prediction error cannot be eliminated to zero.
    • The system is stable (see problem 11.10[1]).
  • Covariance Method only considers n in smaller interval (p to M).
    • There is no end effects like autocorrelation method and thus no need for Hamming window.
    • The prediction error can be zero for deterministic signal. 
    • Stability is not guaranteed.

Professor Oppenheim provided an example in 11.5.3. He tried to estimate the frequency of the sum of two sinusoids in only 101 samples. 

     s[n] = 20cos(0.2*pi*n - 0.1*pi) + 22cos(0.22*pi*n+0.9*pi).

Here I tried to provide a potential matlab code to generate Fig. 11.16 in the DSP [1]. Please note that conventional methods like DFT or PSD cannot distinguish the two components. Neither does the Autocorrelation Method due to the windowing effect. (The main lobe of window is 0.04pi(rad/s) but the frequency difference of the two sinusoids is 0.02pi (rad/s).) Only Covariance Method can separate the two components in this example.



% Code
n= 0:101;
s = 20*cos(0.2*pi*n-0.1*pi) + 22* cos(0.22*pi*n+0.9*pi);

[d, p] = aryule(s,4);
[a, e] = arcov(s,4);

[Ha, wa] = freqz(sqrt(p),d);
[Hc, wc] = freqz(sqrt(e),a);

periodogram(s)
hold on;
hp = plot(wa/pi, 20*log10(2*abs(Ha)/(2*pi)),'r'); % Scale to make on-sided PSD
hp.LineWidth=2;
hp = plot(wc/pi, 20*log10(2*abs(Hc)/(2*pi)),'b'); % Scale to make on-sided PSD
hp.LineWidth=2;
xlabel('Normalized frequency (\times \pi rad/sample)');
ylabel('One-sided PSD (dB/rad/sample)');
legend('PSD estimate of x', 'PSD of model output (auto)', 'PSD of model output (cov)');



[1] Alan V. Oppenheim and Ronald W. Schafer. 2009. Discrete-Time Signal Processing (3rd ed.). Prentice Hall Press, Upper Saddle River, NJ, USA.



Comments

Post a Comment

Popular posts from this blog

[Ch4 ADC/DAC] How to simulate ADC/DAC process in Matlab? How to actually reconstruct a signal nearly sampled in Nyquist Rate?

Image
Reconstruction is essentially a kind of interpolation or so-called digital to analog conversion (DAC). Detail descriptions are introduced in chapter 4.8.3 of the DSP Bible [1]. Although we all understand the basic sampling theory, the problem is how to bridge the theory and practice . Because all the signals we can simulate in Matlab are essential 'Digital'. So, the so-cal ADC/DAC is just an approaching simulation. Let me briefly justify the difference. Ideal Reconstruction in DAC:  There is typically a zero order hold circuit and reconstruction filter. Please note that the frequency response of zero-order hold is not evenly flat so a reconstruction filter comes into an aid. Matlab Simulation in DAC: There is no need to actually simulate a zero-order hold. All you need are an upsampling and an ideal low pass filter. This sample code sample a simulated 2 Hz analog signal,x_a, (with sampling rate 1500 Hz) to a discrete-time signal,x_d, with 5 Hz sampling rate (near
Read more

[Ch7. Filter Design] Why is it a bad idea to filter by zeroing out FFT bins? (為何不能DFT轉過去，直接設成想要的形狀，再IDFT轉回來?)

Image
The complexity of DFT/IDFT using FFT/IFFT is only n*log2n whereas the complexity of convolution is n square. Then, why do we need to learn so many methods to design a filter? Why is it a bad idea to filter by zeroing out FFT bins? This question might be the most confusing one when learning DSP. To verify this fact is simple, we just need to  cautiously  observe the spectrum of the impulse response of an ideal(?) bandpass filter which just zeros out FFT bins. Why do I need to add the adverb "cautiously"?  If we just use the same size of the FFT to observe the response of the impulse, we will be deceived as shown in Fig 1. Nonetheless, if we add the order of DFT when observing the output of the filter, that is, zero-padding the impulse response, we can find the so-called Gibbs phenomenon, ripples in the frequency domain, as depicted in Fig.2. The results in fact come from the windowing effect. If you want to entirely understand the problem, please refer to chapter 7.6
Read more

[Chapter 5.1] Meaning of General Linear Phase or Group Delay of a filter?

Image
For me, one better description of the  general linear phase (GLP)  is a  constant group delay . By definition, group delay is a negative derivative of phase [1] (In fact, the detailed phase/delay relationship can be derived but I suggest we just accept them.). Then, the Derivative (group delay) of a linear (phase) is constant and vise versa. Then, let's move on to the meaning of group delay. the group delay of a frequency represents the delay unit of the filter to that frequency. So, the filter may treat different frequencies with different delay units. For an extremely bad example of a non-linear filter, the input signal 'do re mi' may becomes 're mi do' in the output. A GLP filter can guarantee that such weird condition will never happen. Github link ----------------------------------------------------------------------------------- Here I wrote an example. The example comes from Chapter 5.1.2 in the bible of DSP[2] and I just implemented the sample. Fi
Read more