DSPRelated.com
Forums

Dealing with aliased sinc function when computing spectra

Started by Crandles 2 years ago11 replieslatest reply 2 years ago307 views

Hi,

I currently working on plotting some different spectra of Hamming windows and began exploring their Fourier Transform more. It involves the Aliased Sinc Function or Dirichlet function depending on your nomenclature. This function has a period of 2*pi when it has an odd order and a period of 4*pi when it has an even order. In our case, the window length corresponds to the order. The Discrete-Fourier transform should be periodic of 2*pi, however we often need to use windows of an even length due to the FFT. Taking a closer look at Matlab's functions for generating these windows I noticed there's the option to choose 'periodic' or 'symmetric' when generating the windows, the former attempting to mitigate the issue by computing L+1 samples for the window but only returning the first L. Presumably this is so the signal used for the window actually comes from the spectra of period 2*pi.

It makes sense when you think about it. The L samples returned are going to be from a function which has a spectra which is periodic of 2*pi, however it seems like some of the symmetry in the window is lost. It seems like, depending on you needs in the time-domain this one sample offset and lack of symmetry might be a problem. How much of a difference can it make really though? It mainly seems like if would be an issue if you were modulating the signal close to the Nyquist frequency. This is where I ran into it when I was computing the spectra for an a filterbank interpretation of the FFT.

Any thoughts or experiences people have with the issue would be greatly appreciated.

Thanks,

Graham

[ - ]
Reply by LESheffieldNovember 10, 2021

I personally use an index offset of 1/2.

For example, a Hann of length N (even or odd), normalized for unity point spread function, is given by:

h[n] = 1 - cos(2 pi (n + 1/2) / N)

for 0 <= n < N

I also post-FFT equalize the phase so that it's equivalent to summing over:

(1 - N) / 2 to (N - 1) / 2

instead of 0 to N-1, which is customary.

[ - ]
Reply by Rick LyonsNovember 10, 2021

Hi. Window sequences where the last sample equals the first sample (symmetric) are used in digital filter design. Window sequences where the last sample equals the second sample (periodic) are used in spectrum analysis. Of course the so called "periodic" window sequences are not actually periodic because by definition only infinite-length sequences can be periodic.

[ - ]
Reply by fharrisNovember 10, 2021

Rick,

that's why we examine the signal and the spectrum on the circle where they are periodic.


fred

[ - ]
Reply by Rick LyonsNovember 10, 2021

Hi fred.

Yes, I agree. A finite-length time-domain sequence and its finite-length forward DFT sequence have a special "circularly-periodic" relationship with each other.

Your 8-point window sequence "on a circle" is a fun concept to think about. (Please bear with me.) That sequence can be represented as x(0), x(1), x(2), x(3), x(4), x(5), x(6), and x(7) where x(0) is at the zero-degree point on the circle. If we start at the zero-degree point and go around the circle CCW by 45 degree steps n=7 times (a total angle of 315 degrees) we'll encounter the x(n=7) sample.

Now if we start at the zero-degree point and travel around the circle CCW by 45 degree steps n=9 times (a total angle of 405 degrees) we don't encounter an x(n=9) sample we encounter the x(n=1) sample, or do we? Does an x(n=9) sample exist? For that matter does an angle of 405 degrees exist? I don't have simple "Yes" or "No" answers to those questions.

It seems to me that if you say "Yes, an angle of 405 degrees does exist" then you'd also have to say an x(n=9) sample exists. (I think a shot of Kentucky bourbon, right now, might help me sort this all out.)

[ - ]
Reply by fharrisNovember 10, 2021

Hello Graham,

You never want to use a true even symmetric window. this is a window with a center point of symmetry and a matching left and right sample which gives you an odd number of points. You also want to avoid an even number of points with matching left and right values samples that does not have a center point of symmetry. The problem goes away when you think of the samples on a circle rather than on a line. think of an 8 point window with samples on the circle, 45 degrees apart. You can start at index 0, (0 degrees) at the center point of symmetry and increase the index (+1, +2, +3, +4) to move CCW around the circle as well as decrease the index (-1, -2, -3, -4) to move CW around the circle. Both trips, positive and negative displacement meet at index +4 and -4. But +4 and -4 are the same index and you can't have the same index at both sides of the array. We normally think of the +4 as belonging to the positive frequency span and the -4 is in fact + 4, the start of the next trip around the circle. If you want to see the samples on a line rather than on the circle you take the sample value at +4 and split it, 1/2 its value is at +4 and 1/2 at -4... and you get the odd number of points on the line even though you had an even number of points on the circle. When you wrap the line samples back around the circle the +4 and -4 address samples are added to reconstruct the original. The channelizer wants a window with an odd number of points and zero extended to the left sample: see attached code (it calls myfrf)

Channelizer_20_20_fred_1.m

myfrf.m


See the on line paper "On the Use of Windows for Harmonic Analysis with The Discrete Fourier Transform" (can send me a request for a copy at fjharris@eng.ucsd.edu)


Matlab makes a mistake in their hann (hanning) window design using a true even number of sample without a center value. You know it is incorrect because the fft of the hann window should only be non zero at the indexes (-1, 0, +1 ) or (7, 0, 1)1 but is non zero at every sample except index 4 (for the 8-poit window). In addition, if the hann window were truly even symmetric the imaginary parts would be zero, and they are not.

hanning(8) 

0.1170
0.4132
0.7500
0.9698
0.9698
0.7500
0.4132
0.1170

fft(hanning(8))

 4.5000 + 0.0000i
-1.6941 - 0.7017i
-0.0764 - 0.0764i
-0.0116 - 0.0281i
 0.0000 + 0.0000i
-0.0116 + 0.0281i
-0.0764 + 0.0764i
-1.6941 + 0.7017i    


should have been

hann=0.5-0.5*cos(2*pi*(0:7)/8);

hann

0   
0.1464   
0.5000   
0.8536   
1.0000   
0.8536   
0.5000
0.1464

fft(hann)

 4.0000 + 0.0000i 
-2.0000 + 0.0000i   
 0.0000 + 0.0000i
 0.0000 - 0.0000i
 0.0000 + 0.0000i
 0.0000 + 0.0000i
 0.0000 + 0.0000i
-2.0000 - 0.0000i


fred h

[ - ]
Reply by Rick LyonsNovember 10, 2021

Hi fred.

I've used Matlab's 'hanning()' function for decades. But in my more recent R2019b version of Matlab only a 'hann()' function is documented in their "Documentation" files. However, both 'hann()' and 'hanning()' functions exist in R2019b's "Signal" Toolbox.

windows_39863.jpg

Both the hann(8,'periodic') and hanning(8,'periodic') match the "hann" sequence in your above Reply.


[ - ]
Reply by neiroberNovember 10, 2021

fred and Rick,

Thanks for these very useful comments.  I have been incorrectly using hanning(N) for DFT windowing, when I should have been using hanning(N,'periodic').  I imagine many Matlab users have been making the same mistake, since hanning(N) gives the default symmetric output.

regards,

Neil

[ - ]
Reply by Rick LyonsNovember 10, 2021
Hi Neil. I know exactly what you mean!

Of course the difference between hanning(N,'periodic') and hanning(N,'symmetric') becomes smaller and smaller as N increases.

windows-error_96433.jpg

Niel, I'm a big fan of chebyshev window functions. I just noticed that MATLAB's 'chebwin(N)' command does not produce "periodic" sequences which are needed for spectrum analysis. That's not good! It reinforces my notion that some MATLAB DSP functions were written by mathematicians and not DSP engineers.

[ - ]
Reply by neiroberNovember 10, 2021

Rick,

I vaguely recall fred harris deprecating Chebyshev windows.  I think it was because the flat sidelobes alias.  By contrast, the Kaiser window has sidelobes that diminish with frequency, so aliasing levels are well below the first sidelobe level.  Comments?

Neil


[ - ]
Reply by Rick LyonsNovember 10, 2021

Hi Neil. The reasons I like Chebyshev windows are described in my 13-year old blog (which fred harris helped me write) at:

https://www.dsprelated.com/showarticle/42.php


[ - ]
Reply by CrandlesNovember 10, 2021

Thanks for the replies Fred and Rick. I actually have that paper from a class I took! I had forgotten that we sample the spectrum on the circle and was thinking about things in a very linear fashion. This helps reconcile quite a few other things.