why The FFT can not compute this second harmonic correctly?

[ali_asadzadeh]

Hi,
I have generated two sin waves with the F and 2*F frequency, Suppose that I have a 50Hz signal and a 100Hz signal with equal amplitudes and 64 samples in one cycle of the 50Hz signal

Here is the Matlab code
t=linspace(0,2*pi,64);
x=sin(t);
y=sin(2*t);
sum= x + y;

you can see all the signals at once


Now when I get the full wave FFT with this Matlab command, the results seems ok,

abs(fft(sum))/64

The Result is this,As you can see the first and second harmonic has almost the same amplitude,


Now if I try to cut 32 samples of data (from at some point till 32 samples next) and Calculate only the 32 samples for FFT, I should be able to see the 100Hz amplitude data correctly,since there is at least one cycle of the 100Hz wave data in that window, but As you can see it seems that there is no 100Hz data inside the 32 sample data.

Here is the Matlab code

s= sum(16:47);

And the FFT looks like this

abs(fft(s))/32



what am I missing? Do you have any idea?

[Marco]

You shouldn't use any old FFT to do a DFT, FFT's are often complex FFT ... and to implement DFT you need a "real" FFT (ie. real time domain input -> quadrature frequency domain output). What you're getting out at the moment is complete bullshit. Even if you were doing it right (ie. use DFT instead) you wouldn't see the lowest frequency signal any more in the second case though.

First lets start thinking and expressing the problem more clearly. Don't say FFT, say DFT. Don't say 50 or 100Hz, say frequency is 1/64 and 1/32. When you take 32 samples, the lowest frequency the DFT has a bin for other than DC is 1/32 ... 1/64 simply doesn't have a bin any more.

It's useful to look at the DFT as periodic, you can repeat the block you are transforming and it produces the same DFT result. So what happens when you repeat your 32 sample window? You get lots of sudden jumps at the boundaries, which means high frequency components. All that energy from the low frequency 1/64 has been transferred to high frequencies and DC. Only 1/32 remains, in the first frequency bin.

[awallin]

try again with a time-series that contains 100 or so waves of the thing you want to measure and you will get more sane results.
you'll need an X-axis in Hz to interpret what fft-bin corresponds to what frequency.

numpy/matlab might also have more user-friendly high-level functions for PSD (e.g. "pwelch" or similar)

[rstofer]

Does this help?

http://icrf.nl/Portals/0/Matlab/DFT_and_Matlab.pdf

Lots of math and examples

[dmendesf]

This is a common misunderstanding about what the DFT / DFT does. Most people understand it as "bin N has the exact amplitude / phase of order N harmonic, always". This is just not true. It works like that: "IF the input is a sum of N harmonics that are exactly multiple of the FFT lenght, THEN then bin N has the exact amplitude / phase of the order N harmonic". If you violate this rule then the FFT can have any arbitrary data. I'm posting from my cell phone but when I'm at the computer I'll post a link to a site that has an interactive FFT / DFT Applet that let you change the input vector in real time and look at the output to see what impact it has.

[nick_d]

Firstly, take advantage of superposition i.e. fft(x + y) = fft(x) + fft(y). So do them separately to see what's what. Your 100 Hz that has a full cycle after truncation will show as a spike. Your 50 Hz that has a "hump" after truncation will show as a smear of all frequencies.The hump consists of DC plus a low frequency sine wave phase shifted to have a positive peak in the centre of the hump and negative peaks at the edges. This will give a bit the wrong shape at the edges which it will "correct" by applying more higher frequencies, I think all. You may also find that negative frequencies are important. I don't agree with earlier posting drawing a distinction between DFT and FFT. An FFT is a calculation method for a DFT that gets n log n complexity by a clever grouping of terms. For this kind of work a DFT is just as good, use FFT for performance. The alternative to a DFT is a continuous FT which is only possible using a computer algebra package and does not make sense to apply to vectors of discrete values as here.
cheers, Nick

[ali_asadzadeh]

You shouldn't use any old FFT to do a DFT, FFT's are often complex FFT ... and to implement DFT you need a "real" FFT (ie. real time domain input -> quadrature frequency domain output). What you're getting out at the moment is complete bullshit. Even if you were doing it right (ie. use DFT instead) you wouldn't see the lowest frequency signal any more in the second case though.

First lets start thinking and expressing the problem more clearly. Don't say FFT, say DFT. Don't say 50 or 100Hz, say frequency is 1/64 and 1/32. When you take 32 samples, the lowest frequency the DFT has a bin for other than DC is 1/32 ... 1/64 simply doesn't have a bin any more.

It's useful to look at the DFT as periodic, you can repeat the block you are transforming and it produces the same DFT result. So what happens when you repeat your 32 sample window? You get lots of sudden jumps at the boundaries, which means high frequency components. All that energy from the low frequency 1/64 has been transferred to high frequencies and DC. Only 1/32 remains, in the first frequency bin.

FFT is just the DFT, but we have done it in a fast way

Thanks for the response, so, do you have any idea how I can extract that 100Hz from the 32 sub-samples?

[Dave]

The problem is with using linspace. The spacing between the points isn't 1/64 of 2*pi, but 1/63 of 2*pi. The function spaces the samples so the first and last sample fall directly on the limits that you specified. The frequency of sin(x) therefore is not 1, but 64/63.

Replace the first line with
linspace(0, 2*pi*(63/64), 64)
The FFT is going to show both frequency components being equal.

[ali_asadzadeh]

The problem is with using linspace. The spacing between the points isn't 1/64 of 2*pi, but 1/63 of 2*pi. The function spaces the samples so the first and last sample fall directly on the limits that you specified. The frequency of sin(x) therefore is not 1, but 64/63.

Replace the first line with
linspace(0, 2*pi*(63/64), 64)
The FFT is going to show both frequency components being equal.

Dave thanks, it fixed the first FFT problem and generated a nice 0.5 for 50Hz and 100Hz signals,



But the second half has the energy split in every bin!


and the 100Hz signal has not the peak out there, is there a way to get around this issue?

[Dave]

Half a period of an offset sinewave is not orthogonal to another higher harmonic sinewave, so non-zero results are exactly what you'd expect to see. Try doing the FFT on the same clipped data just with only sin(2x) as the input, you will see just one peak in the 2nd bin (what you, for some reason, call 100Hz).

[GeoffreyF]

A single cycle as you show also has a sharp cut off which generates harmonics.  The ideal FFT is for something that begins infinitely in the past and ends infinitely in the future.  Often a sample signal has many cycles beginning with a very small amplitude and ending the same.

A single cycle is just that, it's not even an approximation of a wave over time.

[Conrad Hoffman]

And thus, windowing, which can be done in about a dozen different ways.

[ali_asadzadeh]

And thus, windowing, which can be done in about a dozen different ways.

Thanks, I was thinking of windowing too, any ideas what type and how should I apply that?

[Marco]

Thanks for the response, so, do you have any idea how I can extract that 100Hz from the 32 sub-samples?

Zero pad the 32 samples to 64, then do the DFT. But you're not really extracting the 1/64 frequency signal, you are assuming it's present and based on that assumption determining it's strength. It could just as well not be, the 32 sample signal is a perfectly valid periodic signal without any 1/64 frequency component at all.

[radioactive]

The problem is that even though the cut samples of 32 length still contains the shorter wavelength, it is still a composite signal.  As others mentioned, the cut causes the boundaries of the finite sequence to be non-zero.  This causes lots of undesired output in the fft.  Windowing can be used to force the boundaries to zero.  I changed the start point of the cut to be 1/8 of the way into the 32-sample sequence and applied triangular windows.  Somewhat of trick question, so here is somewhat of a trick answer:   If nothing else, it shows how to do windowing and zero padding.

Code: [Select]

pkg load signal
N=4096; %zero pad results to N-bin fft
t=linspace(0,2*pi,64);
x=sin(t);
y=sin(2*t);
s1= x + y;
s2= s1((64/8):(64/8)+length(s1)/2-1); %grab half(32) the samples from s1 starting 1/8 in
win1=triang(length(s1)); %generate triangular windows
win2=triang(length(s2));
k1=24  %gain adj for win1
k2=k1/2  %gain adj for win2
s1 = s1.*win1'; %apply windows to samples
s2 = s2.*win2';
s1 = [zeros(1,N/2-32) s1 zeros(1,N/2-32)];  %zero pad for higher-res fft
s2 = [zeros(1,N/2-16) s2 zeros(1,N/2-16)];
fft1 = abs(fft(s1)); %get the mag of fft
fft2 = abs(fft(s2));
plot( fft1(1:length(fft1)/2)/k1, '*r' );
hold on;
plot( fft2(1:length(fft2)/2)/k2, '+b' );
hold off;

[ali_asadzadeh]

The problem is that even though the cut samples of 32 length still contains the shorter wavelength, it is still a composite signal.  As others mentioned, the cut causes the boundaries of the finite sequence to be non-zero.  This causes lots of undesired output in the fft.  Windowing can be used to force the boundaries to zero.  I changed the start point of the cut to be 1/8 of the way into the 32-sample sequence and applied triangular windows.  Somewhat of trick question, so here is somewhat of a trick answer:   If nothing else, it shows how to do windowing and zero padding.

Thanks, Is there away to do it in an arbitrary starting point like the 16:47 points?

[ali_asadzadeh]

Also there is some problems with your solution, Because I'm doing the FFT finally in a Cortex M device, with limited resources and limited time, I have reduced the original signal to 32 samples,so the FFT would be faster and needs less ram resources! Doing it in a 4096 sample windows size would be an over kill, because I have under 1ms to do this stuff on several simultaneous channels of data, so any New Ideas?

[nick_d]

Yes, there are many ways to separate frequencies without FFT. An IIR filter would be an efficient choice. So would a wavelet transform if you need to process a spectrum (but this is usually unnecessary). May we have details of the application?
cheers, Nick

[rstofer]

Here's some MATLAB code I grabbed from the MathWorks site and changed enough to add a couple of harmonics.

Code: [Select]

Fs=1000;
T=1/Fs;
NFFT=256;
FFTperiod=256*T;
t=(0:NFFT-1)*T;
str=['FFT period is ',num2str(FFTperiod), ' seconds'];
Freq=10*1/FFTperiod
strFreq =['Frequency is ', num2str(Freq), ' Hz'];
Freqbin =1/FFTperiod
strFreqbin =['FFT frequency bin is ',num2str(1/FFTperiod), ' Hz'];
strFreqNyquist =['FNyquist is ',num2str(1/FFTperiod*256/2), ' Hz'];
disp(str);
disp(strFreq);
disp(strFreqbin);
disp(strFreqNyquist);
x=sin(2*pi*Freq*t)+(sin(6*pi*Freq*t)/3)+(sin(10*pi*Freq*t)/5);
subplot(2,1,1)
plot(t,x);
title('time domain signal')
Y=fft(x,NFFT)/NFFT;
f=Fs/2*linspace(0,1,NFFT/2);
subplot(2,1,2);
plot(f,20*log10(2*abs(Y(1:NFFT/2))))
title('FFT magnitude in dB20')


[ali_asadzadeh]

Thanks rstofer, Because the application should be so fast, I could only do the FFT in 1 full sin cycle, and for the second harmonic also I wanted to be able to calculate it in half the cycle of the full first harmonic cycle.

[rstofer]

Thanks rstofer, Because the application should be so fast, I could only do the FFT in 1 full sin cycle, and for the second harmonic also I wanted to be able to calculate it in half the cycle of the full first harmonic cycle.

The execution times seem highly variable but the first execution is very slow.

1) 125 ms
2) 798 us
3) 257 us
4) 153 us
5) 512 us

This is on my Surface Book with an I7-6600 2.6 GHz processor.  I imagine it would run a lot faster on my I7-7700 desktop but I have only licensed MATLAB for one machine.

I don't think you are going to get away with just a single cycle of the 2d signal because it needs to be padded with 'something' to become a continuous signal at some integer multiple of the 1st signal.  Any padding other than another cycle of the 2d signal is going to result in harmonics all over the place (I think).

What kind of signal runs as a sin() for 1 cycle and 0v for 1 cycle and then repeats forever?  A composite of a square wave and a sin() wave would do that but we already know that the square wave has odd harmonics from DC to daylight.

[Conrad Hoffman]

There's some good public FFT code available. If you want it fast, code it yourself. http://www.fftw.org/

[ehughes]

Don't forget that but "not windowing"  you are still "windowing" with a rectangular window.        You are seeing the the interaction of the sin(x)/x functions  shifted to the different bins and summed up.

In the case you got all the zeros for most of the bins.   The superposition of the sin(x) is still there,  the bins just line up with nulls in the window function.

In real world cases,   your tones will never line up such that the nulls of the window function are at the bins in your FFT.    You will always be looking at the effects of your window (rectangular or otherwise).  It is just something you have to be aware of when doing your processing.

[Rerouter]

The execution times seem highly variable but the first execution is very slow.

1) 125 ms
2) 798 us
3) 257 us
4) 153 us
5) 512 us

That screams of just in time compilation to me, first run compiles the bytecode, the later only needs to run it.

[ali_asadzadeh]

The execution times seem highly variable but the first execution is very slow.

1) 125 ms
2) 798 us
3) 257 us
4) 153 us
5) 512 us

Thanks, I would run the FFT on an 200MHz cortex M4 part using ARM DSP libraries and I have 8 channels of DATA, so I think I can only pad data to the 128 samples at most, so what window should I choose?

Please note that my ADC is sampling all 8 channels in the ISR with the help of DMA, so I should update each cycle measurements and make some decisions and do some actions based on the results. so regarding all other-things including sampling and the time that's available to me I should find a fast, memory efficient way to do it.