Sinusoidal frequency fixing

Question

This question is based on a previous similar question.

I have the following equation and adjusted (some random data): 0.44 * sin (N * 2 * PI / 30)

I am trying to use FFT to get the frequency from the generated data. However, the frequency ends close, but not equal to the frequency (which makes the wave a little larger than expected)

The frequencies that are maximum for the FFT are 7 Hz, but the expected frequency is (30 / 2PI) 4.77 Hz.

I included a plot of the FFT and plotted values.

The code I use is:

[sampleFFTValues sFreq] = positiveFFT(sampledata, 1); sampleFFTValues = abs(sampleFFTValues); [v sFFTV]= max(sampleFFTValues) 

Positive FFT can be found here . It basically centers the FFT graph and turns off negative signals.

My question is: how can I get the FFT more accurate without resorting to the least squares just for frequency?

Answer 1

I don't think FFT is good for measuring high-resolution frequencies for (quasi) periodic signals - see below.

Each discrete FFT extends to non-integer numbers of bins (that is, at any frequency that does not exactly correspond to one of the frequency stages of a particular FFT); these "intermediate" frequencies will be spread / scattered around the nearest integer bin. The shape of this extension (“extension function”) depends on the window function used for the FFT. This extension feature - to simplify and summarize things - is either very narrow, but very very ragged (very high peaks / very low valleys), or wider but less ragged. In theory, you could do a very thin frequency sweep of sine waves and calculate the FFT for each of them, and then you can “calibrate” the shape and behavior of the function, saving the outputs of all the FFT along with the frequency that led to this conclusion, and then by comparing output FFT signal to be measured, with previously saved results and search for "nearest", you can find a more accurate frequency.

A lot of effort.

But do not do this if you only need to measure the frequency of one signal.

Instead, try measuring the wavelength. It can be as simple as measuring the distance between the zero crosses (perhaps for several cycles to get better accuracy - hell, 1000 cycles if you have a lot of them) in the samples, and divide the sampling frequency by reaching the frequency . Much easier, faster and more accurate.

Example: a sampling frequency of 48000 Hz, a signal of 4.77 Hz leads to a resolution of ~ 0.0005 Hz, simply by measuring the cycle length alone with the roughest approach. (If you execute n cycles, the frequency resolution is also multiplied by n.)

Answer 2

As others have already mentioned, you are misinterpreting the frequency of the signal. Let me give an example to clarify a few things:

 Fs = 200; %# sampling rate t = 0:1/Fs:1-1/Fs; %# time vector of 1 second f = 6; %# frequency of signal x = 0.44*sin(2*pi*f*t); %# sine wave N = length(x); %# length of signal nfft = N; %# n-point DFT, by default nfft=length(x) %# (Note: it is faster if nfft is a power of 2) X = abs(fft(x,nfft)).^2 / nfft; %# square of the magnitude of FFT cutOff = ceil((nfft+1)/2); %# nyquist frequency X = X(1:cutOff); %# FFT is symmetric, take first half X(2:end -1) = 2 * X(2:end -1); %# compensate for the energy of the other half fr = (0:cutOff-1)*Fs/nfft; %# frequency vector subplot(211), plot(t, x) title('Signal (Time Domain)') xlabel('Time (sec)'), ylabel('Amplitude') subplot(212), stem(fr, X) title('Power Spectrum (Frequency Domain)') xlabel('Frequency (Hz)'), ylabel('Power') 

Now you can see that the peak in the FFT corresponds to the original signal frequency at 6 Hz

 [v idx] = max(X); fr(idx) ans = 6 

We can even verify that Parseval's theorem :

 ( sum(x.^2) - sum(X) )/nfft < 1e-6 

---

Option 2

As an alternative, we can use the functions of the signal processing panel:

 %# estimate the power spectral density (PSD) using the periodogram h = spectrum.periodogram; hopts = psdopts(h); set(hopts, 'Fs',Fs, 'NFFT',nfft, 'SpectrumType','onesided') hpsd = psd(h, x, hopts); figure, plot(hpsd) Pxx = hpsd.Data; fr = hpsd.Frequencies; [v idx]= max(Pxx) fr(idx) avgpower(hpsd) 

Note that this function uses a logarithmic scale: plot(fr,10*log10(Pxx)) instead of plot(fr,Pxx)

Answer 3

Assuming N is the time in seconds, your frequency is 1/30 Hz ( y=A * sin( 2* PI * f * t) )

Frequency Resolution = Sampling Rate / FFT Point

The sampling frequency is determined by the nyquist criterion, the sampling frequency (samples / second) should be at least two times the maximum frequency to be analyzed, for example. 48 kHz for analysis up to 24 kHz. (It’s useful to have a bit of buffer for “real life” data).

So, you may need to increase the size of your FFT.

Answer 4

What you are looking for is a frequency estimation method, and there are many. FFT is one component of several assessment methods. Just using a bipolar bean, as in your example, you get the worst resolution (but the greatest noise immunity to any other exactly periodic sinewave). In low noise situations, you can interpolate. Parabolic interpolation of the logarithm value is one common estimate, but synchronization of the interpolation of the FFT results may be better for a rectangular window. Zero filling and performing a longer FFT is basically equivalent to interpolation.

For an exact sinusoid with zero noise, forget about the FFT and just solve the equation in 3 unknowns, which may include only 3 or 4 smoothed sample points, the algorithms for this are here and here .

I list a few other frequency estimation methods on my DSP page .

Answer 5

If you generate a function due to working with samples, you can create many points and run BIG fft, so the frequency bins are very small for high accuracy. But this will not solve the main problem.

Answer 6

Firstly, fixing your question: (30 / 2PI) is not a frequency. Your signal frequency is 1/30 * regardless of the sampling frequency you used. Secondly, can you tell me what is the length of the sampledata vector? When the FFT returns a vector of values, the ith value will correspond to f_i = i / N, where N is the length of the vector and i \ in [0, N-1] You want i / N to exactly equal 1/30 for some integer i . In other words, N should be 30 * i, i.e. N must be a multiple of 30. Now, was the length of the vector you are using a multiple of 30? If you do not try to do this, and that should solve the problem.

Answer 7

Try the window function ?