Continue to Site

Welcome to EDAboard.com

Welcome to our site! EDAboard.com is an international Electronics Discussion Forum focused on EDA software, circuits, schematics, books, theory, papers, asic, pld, 8051, DSP, Network, RF, Analog Design, PCB, Service Manuals... and a whole lot more! To participate you need to register. Registration is free. Click here to register now.

Register Log in

[SOLVED] FFT algorithm and caculate frequency

Status
Not open for further replies.
H

hluyen

Newbie level 5
Newbie level 5
Joined
Mar 13, 2013
Messages
9
Helped
1
Reputation
2
Reaction score
1
Trophy points
1,283
Visit site
Activity points
1,351
Do Anyone known caculate frequency from FFT algorithm with samples input ?
 

For example, 1024pt fft with 44100 sampling frequency using PC.
After doing FFT and found maximum index b put it in function like that:
function Bin2hz(b:integer):integer;
begin
result:=round(b*44100/1024);
end;
 

I previously wrote a script to report the dominant frequency present in a series of captured ADC samples that might suit your need...

It [laboriously] implements a discrete Fourier transform and returns the frequency of the transform bin with the largest power. Note that it was intended to debug an embedded ADC/CPU - it's not written for speed!
Enjoy :)

The (working) MATLAB implementation:

Code C - [expand]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

% Fabricate a test signal to give me something to look at
 
Fs = 1.06E+3;   % Sampling frequency
N = 3125;       % Number of samples
sigF = 189;     % Frequency of test tone
 
sig = sin(2*pi*sigF*(1:N)/Fs) + 0.1*rand(1, N);  % Test tone + noise
 
% Calculate DFT
 
Fr = zeros(1, N);
Fi = zeros(1, N);
 
for k=1:N
    for n=1:N
        Fr(k) = Fr(k)+sig(n)*cos(2*pi*k*n/N);
        Fi(k) = Fi(k)-sig(n)*sin(2*pi*k*n/N);
    end
end
 
% Find peak freq
P = Fr.^2+Fi.^2;
 
peakVal = 0;
for i=1:floor(N/2)
    if P(i) > peakVal
        peakVal = P(i);
        peakLoc = i;
    end
end
 
peakFreq = peakLoc*Fs/N
plot(P);



And a never tested/built C-psuedo code version...

Code C - [expand]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

double primaryFreqComponent(int *ADCdata, int Fs, int N)
  {
  int n, k;
  int peakLoc;            // Location (=frequency) of peak signal component
  double peakVal = 0.0;   // Magnitude of peak signal found
  double Fr[N];           // <--- Buffers of intermediate data the same length as the input data
  double Fi[N];           // (malloc() or statically allocate ... or whatever 
  
  for (k=0; k<N; k++)     // Calculate DFT
     for (n=0; n<N; n++) {
         Fr[k] += (double)(ADCdata[n])*cos(2*pi*k*n/N);
         Fi[k] -= (double)(ADCdata[n])*sin(2*pi*k*n/N);
         }         
 
  for (n=0; n<N; n++)     // Search for peak
     if (Fr(n)*Fr(n)+Fi(n)*Fi(n) > peakVal) {
        peakVal = Fr(n)*Fr(n)+Fi(n)*Fi(n);
        peakLoc = n;
        }
 
  return (double)(peakLoc*Fs)/(double)(N);
  }
 

thylacine1975 said:
I previously wrote a script to report the dominant frequency present in a series of captured ADC samples that might suit your need...

It [laboriously] implements a discrete Fourier transform and returns the frequency of the transform bin with the largest power. Note that it was intended to debug an embedded ADC/CPU - it's not written for speed!
Enjoy :)

The (working) MATLAB implementation:

Code C - [expand]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

% Fabricate a test signal to give me something to look at
 
Fs = 1.06E+3;   % Sampling frequency
N = 3125;       % Number of samples
sigF = 189;     % Frequency of test tone
 
sig = sin(2*pi*sigF*(1:N)/Fs) + 0.1*rand(1, N);  % Test tone + noise
 
% Calculate DFT
 
Fr = zeros(1, N);
Fi = zeros(1, N);
 
for k=1:N
    for n=1:N
        Fr(k) = Fr(k)+sig(n)*cos(2*pi*k*n/N);
        Fi(k) = Fi(k)-sig(n)*sin(2*pi*k*n/N);
    end
end
 
% Find peak freq
P = Fr.^2+Fi.^2;
 
peakVal = 0;
for i=1:floor(N/2)
    if P(i) > peakVal
        peakVal = P(i);
        peakLoc = i;
    end
end
 
peakFreq = peakLoc*Fs/N
plot(P);



And a never tested/built C-psuedo code version...

Code C - [expand]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

double primaryFreqComponent(int *ADCdata, int Fs, int N)
  {
  int n, k;
  int peakLoc;            // Location (=frequency) of peak signal component
  double peakVal = 0.0;   // Magnitude of peak signal found
  double Fr[N];           // <--- Buffers of intermediate data the same length as the input data
  double Fi[N];           // (malloc() or statically allocate ... or whatever 
  
  for (k=0; k<N; k++)     // Calculate DFT
     for (n=0; n<N; n++) {
         Fr[k] += (double)(ADCdata[n])*cos(2*pi*k*n/N);
         Fi[k] -= (double)(ADCdata[n])*sin(2*pi*k*n/N);
         }         
 
  for (n=0; n<N; n++)     // Search for peak
     if (Fr(n)*Fr(n)+Fi(n)*Fi(n) > peakVal) {
        peakVal = Fr(n)*Fr(n)+Fi(n)*Fi(n);
        peakLoc = n;
        }
 
  return (double)(peakLoc*Fs)/(double)(N);
  }
Click to expand...


Thanks !You can give me a fomula caculate frequency from DFT ?I want know from DFT how to caculate frequency .
 

Result of fft or dft is array of amplitudes. Array elements are numbered by indexes. Convert any index to frequency using formula above. It is linear coefficient
 

Why peakFreq = peakLoc*Fs/N . I dont understand.So frequency peakVal and PeakLoc above is what ? Can Termianator3 and Thylacine 1975 explain for me ?
 

#include <stdio.h>
#include <math.h>
#define pi 4.*atan(1.)

I tested with code C below,I see with Fsample 100 is right with result F = 50 Hz but when replace Fsample 200 ,result F = 153 is wrong.I think this code not right.

#define Fsample 100 //sample frequency
#define sigF 50 //signal frequency
#define sampleN 50 // sample number
double primaryFreqComponent(double *ADCdata, int Fs, int N)
{
int n, k;
int peakLoc; // Location (=frequency) of peak signal component
double peakVal = 0.0; // Magnitude of peak signal found
double Fr[N]; // <--- Buffers of intermediate data the same length as the input data
double Fi[N]; // (malloc() or statically allocate ... or whatever

for (k=0; k<N; k++) // Calculate DFT
{
Fr[k] = 0;
Fi[k] = 0;
for (n=0; n<N; n++) {
Fr[k] += (double)(ADCdata[n])*cos(2*pi*k*n/N);
Fi[k] -= (double)(ADCdata[n])*sin(2*pi*k*n/N);
}
}
for (n=0; n<N; n++) // Search for peak
{
if (Fr[n]*Fr[n]+Fi[n]*Fi[n] > peakVal) {
peakVal = Fr[n]*Fr[n]+Fi[n]*Fi[n];
peakLoc = n;
}
}
return (double)(peakLoc*Fs)/(double)(N);
}

int main()
{
double ADCsample[sampleN];
double F,sampletime;
int i;

for(i = 0;i<sampleN;i++){
sampletime = (double)i/Fsample;
ADCsample = sin(2*pi*sigF*sampletime);
}
F = primaryFreqComponent(ADCsample, Fsample, sampleN);
printf("Freq :%f",F);
}
 
Well, that turned out to be an insidious little bug! The good news is the code is not fundamentally flawed, it just needs a tweak :)

The problem is in the loop that searches for the peak amplitude frequency bin - for 'simplicity', I've searched over every transformed frequency bin (i.e. all n<N). This isn't necessary, since for a real valued input (such as a single stream of ADC samples) the DFT bins above N/2 are (theoretically) the mirror image of the bins below N/2. [The frequency of bin N/2 corresponds to the Nyquist frequency of the specified sampling rate].

As luck (!) would have it - due to numerical rounding errors/precision - for the values you specified, the magnitude of (image) bin #38 was ~1E-9 larger than it's (fundamental) in bin #12. As such, the frequency that was being reported was of the image, i.e. Fs-50 = 150Hz... The 3Hz error (i.e. 153 Hz instead of the correct 150Hz) is a consequence of the coarse frequency binning - with N=50, the frequency resolution is only ~8Hz (~= 2Fs/N).

The fix is easy though! Simply change the loop declaration/conditional to only search the frequency bins below Nyquist. i.e.:

Code C - [expand]
1
2

for (n=0; n<(N/2); n++) // Search for peak
 { ...


Works like a charm now... enjoy!
(...and thanks for adding a wrapper/testing it! Neat :)


Come to think of it - if we're not going to search over the bins > N/2, then there's no point calculating them either. The previous loop involving k could also be similarly truncated:

Code C - [expand]
1
2

for (k=0; k<(N/2); k++) // Calculate DFT
 { ...


(...as well as the declarion of the arrays Fr and Fi...)
 
Last edited:

Peakfrequency is Peakloc x fs / n. For example you use n 1024 point fft. If you observe n 1024 point signal, the shortest sine period is 2 points, one in the bottom and one in the top. So maximum frequency is 512 two-point sines per 1024 points. And peakloc gives your maximum index of 512. Dividing peakloc by n you can obtain percentage of maximum frequency. Next multiply it by your sample frequency. There is no way to obtain sample frequency from signal, it is know from hardware. For computers sampling frequency fs 44100 is common. So the formula is peakloc x 44100 / 1024
 

thylacine1975 said:
Well, that turned out to be an insidious little bug! The good news is the code is not fundamentally flawed, it just needs a tweak :)

The problem is in the loop that searches for the peak amplitude frequency bin - for 'simplicity', I've searched over every transformed frequency bin (i.e. all n<N). This isn't necessary, since for a real valued input (such as a single stream of ADC samples) the DFT bins above N/2 are (theoretically) the mirror image of the bins below N/2. [The frequency of bin N/2 corresponds to the Nyquist frequency of the specified sampling rate].

As luck (!) would have it - due to numerical rounding errors/precision - for the values you specified, the magnitude of (image) bin #38 was ~1E-9 larger than it's (fundamental) in bin #12. As such, the frequency that was being reported was of the image, i.e. Fs-50 = 150Hz... The 3Hz error (i.e. 153 Hz instead of the correct 150Hz) is a consequence of the coarse frequency binning - with N=50, the frequency resolution is only ~8Hz (~= 2Fs/N).

The fix is easy though! Simply change the loop declaration/conditional to only search the frequency bins below Nyquist. i.e.:

Code C - [expand]
1
2

for (n=0; n<(N/2); n++) // Search for peak
 { ...


Works like a charm now... enjoy!
(...and thanks for adding a wrapper/testing it! Neat :)


Come to think of it - if we're not going to search over the bins > N/2, then there's no point calculating them either. The previous loop involving k could also be similarly truncated:

Code C - [expand]
1
2

for (k=0; k<(N/2); k++) // Calculate DFT
 { ...


(...as well as the declarion of the arrays Fr and Fi...)
Click to expand...

I still this code has problem when fixed:
for (n=0; n<(N/2); n++) // Search for peak
{ ...

result after fixed and redefine

#define Fsample 400 //sample frequency
#define sigF 200 //signal frequency
#define sampleN 100 // sample number

Result of Frequency is 24Mhz.

Please,You can review!
 

24 MHz eh? How on earth did it come up with that?!
I don't have my compiler handy at the moment, but I can guess - you've given it a pretty testing set of input parameters... an input frequency *exactly*=Nyquist (SigF = Fsample/2).
Whether or not the peak search even finds this peak will depend on whether sampleN is even or odd, since the loop terminates when n<(N/2).

Try with a slightly friendlier input signal and see how it goes :)

P.S. I just thought of another reason this mightn't be working: With SigF = Fsample/2, your test input vector:
Code C - [expand]
1

ADCsample[i] = sin(2*pi*sigF*sampletime);

...will be zero for all values of i!
 

Sorry Frequncy F is 24Hz.Why with sigF = Fsamle/2,all valules of vector ADCsample is zero?I debuged and is not zero.Theoretically of Nyquist Fsample = sigF*2.With ADCsample is sampled at 200hz frequency above,then F result after DFT must 200hz instead 24hz.
 

Hi all!use FFT can measure square wave frequency ? Or only measured sinwave.?
 

Square wave will give many harmonics and leakage. For square wave simpler method will give better result, something like period histogram.
 

Hi guys!

Aye - as Terminator3 said, there are a handful of other (easy) techniques for frequency measurement you can use. Don't despair if you're only using the DFT approach though - it'll cope with any waveshape (sine/square/triangle/mixture etc) you throw at it and simply report the frequency of the dominant component. To illustrate this, I fed the routines listed above a square wave of a specified frequency...

(which I created by generating a sinusoid in the same fashion as you did hluyen, and then "squaring it up" by:
Code C - [expand]
1
2
3
4
5
6

for(i = 0;i<sampleN;i++){
if (ADCsample[i] > 0)
   ADCsample[i] = 1;
else
   ADCsample[i] = 0;
}


Plotting the magnitude of the computed DFT bins shows the expected odd harmonics of a square wave:


...but most importantly shows that the frequency was correctly estimated. Note that I (randomly) selected some horrible Fs/N/Fsig parameters to ensure there were no "friendly" integer/harmonic relationships between them in order to exaggerate the effects of spectral leakage (which is pretty minor in this example, and won't hurt the accuracy of your frequency estimates anyway [finite frequency bin width is the limiting factor]).

- - - Updated - - -

Oh haohaodk46, I wasn't sure what you were asking of me...? Are you seeking other references? If it's code that you're after, hluyen's post #7 can be cut and paste into an editor and fed straight to a compiler such as gcc with (for example): gcc -lm [filename.c] -o [yourOutputFile]
 

I modify code that get sample of sinwave 50 Hz as below :

for(i = 0;i<sampleN;i++){
temp = 1/((double)sigF*sampleN);
sampletime = (double)i*temp;
ADCsample = sin((double)2*pi*sigF*sampletime);
}

and define Fsample must equal 1/((1/sigF)/N).example with signalF = 50Hz then Fsample = 1/(1/50)/8) = 400 hz.Code alway right when increase sample point number and singF with Fsample as formula above.
 

Hi,
I don't see any word about "windowing" in FFT. In case you don't analyze non-periodic real time signals the formula listed above will be OK. But once you start analyzing real time signals (continues stream of data) you need to apply "windowing" function otherwise your output results will be strongly distorted. The idea of FFT is to transform a piece of periodic signal, this piece of signal should also repeat again and again. Of course FFT could be used over non-periodic signals but only after applying windowing functions. There are many known windowing functions - each has some advantages and disadvantages.
 

Absolutely - windowing is required for any sensible spectral estimation - but in this case the aforementioned DFT routine is solely intended to find the dominant frequency present in a bunch of ADC samples (where windowing is overkill, and spectral leakage doesn't cause frequency estimation errors). The plot of bin magnitudes in post #16 was for illustrative purposes only - the function(s) of posts #3,#7 only return a single value.

I wouldn't recommend using the "manually computed" DFT above for "real" spectral estimation - it's much too slow for starters. There's a stack of optimised FFT routines available (one of my favourites is FFTW: https://www.fftw.org/). These, combined with windowing/overlap (e.g. https://en.wikipedia.org/wiki/Welch's_method) will give you much better results!
 

Hi,
I agree - if there is a strong dominant frequency windowing is not needed. If the dominant frequency strength is close to other frequencies then windowing is required. In all cases it's advisable to reduce the calculations, so if we can avoid the windowing - we should avoid it.
 

Status
Not open for further replies.

Similar threads

Part and Inventory Search

Welcome to EDABoard.com

Sponsor

Back
Top