Normalisation of FFT not consistent.

Hi, I am having a few issues normalising my frequency and amplitude for a test FFT I am trying to code. I think it boils down to my understanding of the relationship between signal length, sampling frequency etc.

Below is an example of my code for one set of configuration parameters. We are looking to convert the frequency space results into the physical space. Hence I want to recover a 20Hz frequency and an amplitude of 3 from a FFT.

clc; clear all; close all;

%%Creating test data to ensure FFT is working correctly
Ls = 10;       % Length of signal
ds = 10000     % discrete data points
var = linspace(0,Ls,ds);

% Signal Frequencies and amplitudes
kx = 20; ax = 3;

% Signal creation
D = ax*cos(2*pi*kx*var);

%%FFT Stage
Fs = 1000;   % Sampling Frequency
n = 2^nextpow2(Fs); 
Y = fft(D,n);

% Normalising Y amplitude by number of sampling points
Y = abs(Y/n);

%%Plotting
% Creating vector f that transforms FFT frequency bins to frequency space
% defined by var.
f = linspace(0,1,length(Y))*Fs;

plot(f,Y)
xlabel('Frequency')
ylabel('Amplitude')

This produces the plot:

Which is giving me the correct frequency but not the right amplitude (I am not concerned with the one sided spectrum at the moment).

If I increase my sampling rate by an order of magnitude to Fs = 10000, my frequency changes and my amplitude is still not correct.

I have tried a few permutations of scaling behaviour, but cannot isolate the exact problem. What is the relationship between my signal length (Ls), the resolution of the data over that length (ds), and the sampling frequency used for FFT of that data (Fs)?

Also, how do I best go about increasing the accuracy of FFT so that my frequency is as close and sharp at 20Hz with magnitude exactly 3?

Please help me signal processing gurus! I am eager to learn!

Answer 1

[Star Strider]

Three problems:

n = 2^nextpow2(ds);                                         % ← Use ‘ds’ Not ‘Fs’
Y = fft(D,n)/ds;                                            % ← Divide by ‘ds’ Here

Y = abs(Y);                                                 % ← Do Not Divide By ‘n’ Here

You have to divide by the length of the original signal, because that is where all the energy is. There is no additional energy in the zero-padded signal.

That should solve the amplitude problem. Note that the signal energy is divided into two equal peaks in the two-sided Fourier transform, so the amplitudes will be half the original signal amplitude. In the one-sided Fourier transform, you need to multiply the amplitude by 2.

[olski one]

Thanks! That clears up a lot of my stupid mistakes.

[Star Strider]

My pleasure.

I wouldn’t call them ‘stupid’. It takes time to understand how the Fourier transform works.

[David Goodmanson]

The process still does not result in sharp, one-point peaks with amplitude 3 (if doubled), though.

[Star Strider]

Of course not. It’s a numeric approximation.

[David Goodmanson]

that's not what I meant. With this code the peak maxes out at about 1.425, and there are contributions down at the bottom of the peak at frequencies other than 20. With a combination of a better initial waveform and not using nextpow2, you can end up with two one-point peaks at 20 and -20, amplitude 1.5, zero everywhere else, all good to many many decimal places.

[Star Strider]

My point exactly.

[David Goodmanson]

really? well, in light of his question:

"Also, how do I best go about increasing the accuracy of FFT so that my frequency is as close and sharp at 20Hz with magnitude exactly 3?"

it seems unlikely that the OP realizes that one could have done better, that's all.

[olski one]

I did realise I could have done better. Which was really at the crux of my question you highlighted.

But thank you for the input David. It helped immensely.

So what I have gathered so far is (and please correct if I am wrong):

• The number of DFFT points (n) should be at least 2*ds (the number of data points) -- For accuracy of the FFT.
• We use the nextpow2 to pad zeros not for accuracy of FFT but for speed of the algorithm
• If we use the n = ds on an input signal which is exactly periodic and coherently sinusoidal, then the peaks will be exactly at the input frequency and half the input magnitude.
• In reality however, we do not know the nature of the input signal, and it is not likely to be exactly periodic and coherent, hence, padding zeros should not detract from isolating the dominant frequency and magnitudes.
• Lastly, if considering FFT in the spatial domain. The frequency is refered to as the wavenumber. Each result in the FFT spectrum is the magnitude of the Fourier coefficient (C_k) at the discrete wavenumber (k).

Does that all sound okay? I have only started to teach myself this topic yesterday so trying to wrap my head around it.

[David Goodmanson]

Hi olski,

I don't think any of these points are incorrect except for probably the second one, and I don't understand the first one. Some of this is just difference of opinion, but I have my own take so here are some comments, which ran long because of [2]. Except for the last point I'll use time and frequency as the parameters.

[1] There is the Nyquist rule that you need at least two sample points per oscillation in the time domain to catch a frequency, which is the absolute idealized minimum. For quality results with real data you need more than that. But what does 'two fft points per data point' signify?

[2] I think that in 2017 this idea is mostly folklore. I'm pretty sure that modern fft routines know how to deal efficiently with factors of 3,5, etc. in the number of points. People are so inured to 2^n though, that it's used incessantly. Many times people will pad a tiny 1000 point array to 1024. I wish I knew why. The rule for an fft is always that for the time and frequency array spacings delt and delf,

delt*delf = 1/N.

This works great for N = 1000 but if N is 1024 then at least one of the array spacings is quite awkward. And the execution time gained, if any, is infinitesimal.

For any N, if the signal is periodic in the original time array, then as you know padding with zeros messes that up and affects the frequency domain, sometimes not too much but sometimes unacceptably.

All right, let's go to a lot more points so 2^n can show its stuff. For larger arrays you don't want to have a prime number of points since that's slow, but consider, say, a million points. 2^20 is about 5% more than 1e6, so 2^20 is at a slight disadvantage, but its increased speed should more than make up for it, right? However, on my pc, for a complex array the 2^20 fft is about 5% slower than the 1e6 fft, and for a real array the 2^20 fft is about 17% slower. Going up to even more points, around 1e8, does not change things. It's not hard to find an N where 2^n has 90% as many points as N but is 15% slower.

Where is the vaunted speed increase?

I think this effect may have to do with the amount of microprocessor cache available, but it's a real world result and I think my pc provides a normal example.

[3] I agree, although your code even with the Star Strider modifications is not an example of that.

[4] In terms of dominant frequencies I agree. For magnitudes you have to be a bit careful since the zero padding affects them. But for the most part I think it works, given that you believe padding with zeros is advantageous for some reason or other.

[5] Yes, although there is this terminology problem with 'wave number'. For time and frequency there is the period T, and f and w. For space you only have wavelength lambda, and wave number k. Then

              T * f   = 1            T * w = 2*pi
         lambda * (?) = 1       lambda * k = 2*pi

That's how k is defined in physics, anyway. There doesn't seem to be a good generic term for inverse wavelength. Spectroscopists use inverse cm but that's unit-dependent.

[olski one]

Thanks for the really in depth reply. What I gathered most was that, like most things, they should always be taken with a grain of salt.

In response to some of the points

[1] I observed that suggestion being used in a few examples I was reading. The way I justified it, and I suppose consistent with the Nyquist criteria, the highest frequency captured by the data is that between two adjacent data points. Hence, if we want to adequately capture that smallest period (per Nyquist), then 4 DFFT points would be needed as a minimum. Hence, 2*the number of data points.

[2] Thanks for clarification on this. What I take is that it is not a hard and fast rule and situation dependent. Caveats, Caveats, Caveats!

[3] Why is my example with Star Strider's modifications not an example of a perfectly sinusoidal signal on a periodic domain?

[David Goodmanson]

[1] I still think this is not resolved. The two points per cycle rule is for the maximum frequency that can (barely) be caught from the data. It's a limitation in what the data can provide. If all the frequency components in the data are low enough that they have several points per cycle, the fft is good as is. If there were a frequency large enough that is has less than two points per cycle, effectively it isn't in the data. The time sampling wasn't quick enough to capture it. Doing some kind of 2x procedure isn't going to get it back. You would have to go back and sample the data at a higher rate.

[2] I was trying to emphasize that 2^n may be situation dependent, but I can't really think of a situation where you would ever want to use it, unless you had a large data set with a prime number of points. A comment on [3] will be a bit later, but for now I'll just mention that nextpow2 had a major role in screwing up your code so it's not periodic.

[David Goodmanson]

I agree that with real data you can usually get away with a less exacting approach as in [4], not that that is something to be proud of. But your original point [3] specified some very precise conditions and here is some demo code for that. I used sine instead of cosine because the flaws are easier to see.

By definition, every fft is periodic in whatever you put into its input array. Whether it's periodic in what you want is another question. You can see what the periodicity looks like by concatenating a few arrays side by side.

Here is your coding method for a smaller size array, concatenated. For the plot, don't worry about the x axis scaling, just the shape. At 50 and 100 you see jogs. That's because with linspace, sin = 0 at both ends of the array, so the zero is repeated side by side. You are not periodic in a true sine wave, and that messes up the acuity of the fft.

N = 50;   T = 1;   f0 = 2;     % T is time duration
t = linspace(0,T,N);
y = sin(2*pi*f0*t);
figure(1)
plot([y y y])

The idea is, at the end of the array stop one point short of sin equaling 0, so that you get a true repeating sine wave per the next plot. If you fft this you will get two points with spikes, zero everywhere else.

N = 50;   T = 1;   f0 = 2;     % T is time duration
t = ((0:N-1)/N)*T;
y = sin(2*pi*f0*t);
figure(2)
plot([y y y])

Here is a reconstruction of what nextpow2 does to the function. If using linspace causes some problems, here the periodicity of the sine is totally whacked.

N = 50;   T = 1;   f0 = 2;     % T is time duration
t = ((0:N-1)/N)*T;
y = sin(2*pi*f0*t);
ypad = zeros(1,64);
ypad(1:50) = y;
figure(3)
plot([ypad ypad ypad])

Your present modified code stops short of peak amplitude 1.5 and has broadened out peaks and not spikes. If you make appropriate changes and exile nexpow2 to the outer darkness you can get two spikes at f = +-20, height = 1.5, all with excellent accuracy.

[olski one]

I really cannot thank you enough for the input David. Your examples make your points very clear. These are all things that are not that easy to pick up in the literature.

I feel like I am much better equipped to begin my journey into Fourier space and beyond!

Are there any textbooks you would recommend for clarity and comprehension?

Answer 2

[David Goodmanson]

Good question. I didn't take a course in DSP, so I don't have a textbook recommendation from experience. Rodriguez mentions the problem about the one extra array point at least twice (although not in terms of a linspace function) but he does not really emphasize it. The book has a lot of good information but it is 30 years old, kind of lab oriented, and pretty dense in the text. I wouldn't recommend rushing off and buying a copy. I have not looked at the Schaum's Outlines DFT book but Schaum's is usually pretty good (their complex variables and topology books are excellent), they have lots of solved problems and the price is right. I'm sure there must be a library you can go graze in, but if I had to buy a book sight unseen that would be it.