{"id":3819,"date":"2020-07-03T09:01:12","date_gmt":"2020-07-03T09:01:12","guid":{"rendered":"http:\/\/www.raymaps.com\/?p=3819"},"modified":"2024-12-09T08:17:37","modified_gmt":"2024-12-09T08:17:37","slug":"a-comparison-of-fft-music-and-esprit-methods-of-frequency-estimation","status":"publish","type":"post","link":"https:\/\/www.raymaps.com\/index.php\/a-comparison-of-fft-music-and-esprit-methods-of-frequency-estimation\/","title":{"rendered":"A Comparison of FFT, MUSIC and ESPRIT Methods of Frequency Estimation"},"content":{"rendered":"\n

As discussed in previous posts it is\nfrequently required in communications and signal processing to estimate the\nfrequency of a signal embedded in noise and interference. The problem becomes\nmore complicated when the number of observations (samples) is quite limited.\nTypically, the resolution in the frequency domain is inversely proportional to\nthe window size in the time domain. Sometimes the signal is composed of\nmultiple sinusoids where the frequency of each needs to be estimated separately.\nSimple techniques such as Zero Crossing Estimator fail in such a scenario.  Even some advanced techniques such as MATLAB\nfunction \u201cpwelch\u201d fail to distinguish closely spaced sinusoids. <\/p>\n\n\n\n\n\n\n\n

In this post we discuss two of the most\npopular subspace methods of frequency estimation namely MUSIC* and ESPRIT**. As\na reference FFT method is also presented. Both MUSIC and ESPRIT work by\nseparating the noisy signal into signal subspace and noise subspace. MUSIC\nworks by exploiting the fact that noise eigen vectors that compose the noise\nsubspace are orthogonal to the signal vectors (or steering vectors). So, we can\nsearch for the signal vectors that are most orthogonal to the noise eigen\nvectors and that gives us a frequency estimate. MUSIC is in fact an advanced\nform of Pisarenko Harmonic Decomposition (PHD) where the model order is one\nmore than the number of sinusoids. There is no such limitation in MUSIC. <\/p>\n\n\n\n

Unlike MUSIC which exploits the noise subspace, ESPRIT method uses the signal subspace. It estimates the signal subspace S from the estimate of the signal correlation matrix R (same correlation matrix that was used to estimate noise subspace in MUSIC algorithm). Here is the trick you need to perform on the eigen vectors forming the signal subspace S. Split the matrix S into two staggered matrices S1 and S2 of size (M-1) x p each, where M is the model order and p is the number of sinusoids. S1 is matrix S without the last row and S2 is the matrix S without the first row. Now divide the second matrix S2 by S1 using the Least Squares (LS) approach to obtain matrix P. The angles of eigen values of P give us the estimate of the angular frequency.  <\/p>\n\n\n\n

It is seen that frequency resolution of FFT and MUSIC methods depends upon the size of the temporal window. Greater the length of the time window higher is the frequency resolution. In the results we have shown the frequency resolution is 4MHz. But there is no such barrier for the ESPRIT method. In fact, in the absence of noise, we have observed that ESPRIT method can even decipher sinusoids placed 0.1 MHz apart. But once noise is added the frequency estimation of ESPRIT is not that great. So, for this method the limiting factor is the AWGN noise added to the signal not the size of the temporal window. Lastly, for the simulation scenario we have considered the model size M for MUSIC and ESPRIT is set to 100. We would further investigate this in a future post. <\/p>\n\n\n\n

\"\"
PSD of Two Tones 20 MHz Apart, SNR=7dB, SIR=10dB<\/figcaption><\/figure>\n\n\n\n
\"\"
PSD of Two Tones 8 MHz Apart, SNR=7dB, SIR=10dB<\/figcaption><\/figure>\n\n\n\n
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n%       SPECTRAL ESTIMATION USING\n%      FFT, MUSIC\/PISARENKO, ESPRIT\n%           www.raymaps.com\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\nclear all\nclose all\n\nN=1001;             %Number of samples\nf1=1.0000e9;        %Frequency of tone 1\nf2=1.0200e9;        %Frequency of tone 2\nfs=4*f1;            %Sampling frequency\nw1=2*pi*f1\/fs;      %Angular frequency of tone 1\nw2=2*pi*f2\/fs;      %Angular frequency of tone 2\n\n%%%%% Signal Generation in AWGN Noise %%%%%\nn=0:N-1;\ns1=sqrt(1.00)*exp(i*w1*n);\ns2=sqrt(0.10)*exp(i*w2*n);\nwn=sqrt(0.10)*(randn(1,N)+i*randn(1,N));\nx=s1+s2+wn;\nx=x(:);\n\n%%%%% FFT Spectrum Generation %%%%%\nf=0:fs\/(N-1):fs;\nFFT_abs=abs(fft(x));\nplot(f,20*log10(FFT_abs\/max(FFT_abs)),'linewidth',3,'b+-');\nhold on\n\n%%%%% MUSIC\/Pisarenko Spectrum Generation %%%%%\nM=100;\np=2;\n\nf=0:fs\/(N-1):fs;\nw=2*pi*f\/fs;\n\ncov_matrix=zeros(M,M);\nfor n=1:N-M+1\n  cov_matrix=cov_matrix+(x(n:n+M-1))*(x(n:n+M-1))';\nend\n\n[V,lambda]=eig(cov_matrix);\ne=exp(j*((0:M-1)')*w);\n\nden=zeros(length(w),1);\nfor n=1:M-p\n  v=(V(:,n));\n  den=den+(abs((e')*v)).^2;\nend\nPSD_MUSIC=1.\/den;\nplot(f,10*log10(PSD_MUSIC\/max(PSD_MUSIC)),'linewidth',3,'ro-')\nhold on\n\n%%%%% ESPRIT Based Frequency Calculation %%%%%\nM=100;\np=2;\n\ncov_matrix=zeros(M,M);\nfor n=1:N-M+1\n  cov_matrix=cov_matrix+(x(n:n+M-1))*(x(n:n+M-1))';\nend\n\n[U,E,V]=svd(cov_matrix);\nS=U(:,1:p);\nS1=S(1:M-1,:); \nS2=S(2:M,:);\nP=S1\\S2;\nw=angle(eig(P));  \nf_est=fs*w\/(2*pi); \n\n%%%%% Plotting the PSD %%%%% \nline([f_est(1), f_est(1)],[-40, 0],'Color','Black','LineStyle','--','linewidth',3);\nline([f_est(2), f_est(2)],[-40, 0],'Color','Black','LineStyle','--','linewidth',3);\ntitle('Spectrum')\nxlabel('Frequency')\nylabel('Power')\nlegend('FFT','MUSIC','ESPRIT')\naxis([0.98e9 1.04e9 -40 10])\ngrid on \nhold off\n<\/pre>\n\n\n\n
Note:<\/strong> <\/p>\n\n\n\n
  1. ESPRIT method can be thought of as differential demodulation of a Continuous Phase Modulated (CPM) signal. But instead of working with staggered samples we are working with staggered matrices. It\u2019s the rotation that want to measure in both.<\/li>
  2.  One point to be noted is that both MUSIC and ESPRIT methods require the knowledge of number of sinusoids embedded in noise. FFT on the other hand does not need this information beforehand. <\/li>
  3. From the code it may seem that noise variance is 0.1 but its actually 0.1 per dimension. The total noise variance is 0.2 resulting in an SNR of 10*log10(1\/0.2)=7dB.<\/li>
  4. At a Sampling Frequency of 4GHz, 1001 samples translate to a time window of only 250nsec.  Thus the resolution of FFT method is 1\/250nsec=4MHz.  <\/li>
  5. Please note that the Zero Crossing algorithm we presented in the last post can be easily modified to work with complex exponentials. After all, a complex exponential can be broken down into a real and imaginary part and Zero Crossings of each can be calculated separately.   <\/li><\/ol>\n\n\n\n
    *MUSIC: Multiple\nSignal Classification<\/p>\n\n\n\n
    **ESPRIT: Estimation\nof Signal Parameters via Rotational Invariant Techniques<\/p>\n\n\n\n
    References:<\/strong><\/p>\n\n\n\n
    https:\/\/en.wikipedia.org\/wiki\/MUSIC_(algorithm)<\/a><\/p>\n\n\n\n
    https:\/\/en.wikipedia.org\/wiki\/Pisarenko_harmonic_decomposition<\/a><\/p>\n\n\n\n
    https:\/\/en.wikipedia.org\/wiki\/Estimation_of_signal_parameters_via_rotational_invariance_techniques<\/a><\/p>\n","protected":false},"excerpt":{"rendered":"
    As discussed in previous posts it is frequently required in communications and signal processing to estimate the frequency of a signal embedded in noise and interference. The problem becomes more complicated when the number of observations (samples) is quite limited. Typically, the resolution in the frequency domain is inversely proportional to the window size in the time domain. <\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[81],"tags":[224,82,70,223,110],"class_list":["post-3819","post","type-post","status-publish","format-standard","hentry","category-fundamentals","tag-esprit","tag-fft","tag-frequency","tag-music","tag-psd"],"_links":{"self":[{"href":"https:\/\/www.raymaps.com\/index.php\/wp-json\/wp\/v2\/posts\/3819","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.raymaps.com\/index.php\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.raymaps.com\/index.php\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.raymaps.com\/index.php\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.raymaps.com\/index.php\/wp-json\/wp\/v2\/comments?post=3819"}],"version-history":[{"count":13,"href":"https:\/\/www.raymaps.com\/index.php\/wp-json\/wp\/v2\/posts\/3819\/revisions"}],"predecessor-version":[{"id":4540,"href":"https:\/\/www.raymaps.com\/index.php\/wp-json\/wp\/v2\/posts\/3819\/revisions\/4540"}],"wp:attachment":[{"href":"https:\/\/www.raymaps.com\/index.php\/wp-json\/wp\/v2\/media?parent=3819"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.raymaps.com\/index.php\/wp-json\/wp\/v2\/categories?post=3819"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.raymaps.com\/index.php\/wp-json\/wp\/v2\/tags?post=3819"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}