FFTW Amplitude and Frequency

[skiner36]

Hi,
I am a C-Sharp developer and a beginner at C++.

I am using C++ to get learn about FFTW but unfortunately I haven’t been able to find any examples of FFTW usage that I understand. The sort of thing that makes the lightbulb go on in my head.

I was hoping that someone here might be able to help me understand how to do something that I believe should be quite simple with FFTW.

I have an one dimensional array of double values and I want to find two lots of information from it. One is the cycle amplitude and the other is the cycle frequency.

I have searched this site but I haven’t been able to find anything that fits. I am not an electronics engineer so a lot of what I read on FFTW.org goes over my head a bit.

I should also say that I am reading the values into the array from a text file one line at a time and adding the contents to the array then I am writing the calculated values to another text file like follows. Input value, Frequency, Amplitude.

writing the
Any help would be greatly appreciated

Regards,

John

[ne555]

std::vector<double> in(n); //your data in time domeain (note, real)
std::vector<std::complex<double> > out(n/2+1); //its fourtier transform (note, complex)


fftw_plan p = fftw_plan_dft_r2c_1d( //this would be like a constructor
	in.size(),
	in.data(),
	reinterpret_cast<fftw_complex*>(out.data()), //fftw3 uses its own data type, need to cast
	FFTW_ESTIMATE
);

//read the data from file
//don't exceed the size, as that would need a reallocation
//according to the docs «You must create the plan before initializing the input»
//may use an temporary instead
/*read(foo)
in.resize(foo.size());
fftw_plan p(...)
in = foo*/

fftw_execute(p); //performs the fourier transform, fills `out'


//operate on the fourier transform
for(int K=1; K<out.size(); ++K){ //starts in 1 to ignore DC
	if(std::abs(out[K]) > threshold)
		std::cout << "Frequency " << K*delta_f << " Amplitude " << 2*std::abs(out[x])/in.size() << '\n';
}

fftw_destroy_plan(p); //destructor 

some things to notice
- fftw3 is a C library so it doesn't know about c++ containers like std::vector or std::complex. Instead it provides its own data types and functions for complex numbers and memory management, that's why the casting is necessary.
- also because of C, needs an explicit call to the destructor
- your use case were real signals in time, that gives you conjugate complex in frequency. Because of that some optimizations can be made and will only give you half of the transform.
- «FFTW computes an unnormalized DFT» that's why you need to divide by the number of samples: in.size(). (the factor 2 is because you are only seeing one half)

[skiner36]

Hi ne555,

Thank you very much for your reply it is a big help.

I have a couple of questions if I may:

1. Could you explain what delta_f is and how I would calculate it?
2. Is threshold the minimum value that I want to report on?

Thanks again,

John

[Ganado]

Each bin (index) of the discrete fourier transform is a frequency, but because the DFT just works with a raw array of numbers, you the user have to know what each frequency bin actually means, based on what you know the frequency the samples were sampled at.

See: https://electronics.stackexchange.com/questions/12407/what-is-the-relation-between-fft-length-and-frequency-resolution

I'll just quote it,

The frequency resolution is dependent on the relationship between the FFT length and the sampling rate of the input signal.

If we collect 8192 samples for the FFT then we will have:

                   8192 samples = 4096 FFT bins                       

If our sampling rate is 10 kHz, then the Nyquist-Shannon sampling theorem says that our signal can contain frequency content up to 5 kHz. Then, our frequency bin resolution is:

                   5 kHz/4096 FFT bins ≃ 1.22 Hz/bin                      

The delta_f is the Hz/bin.

Unit-wise: (bin) * (Hz/bin) = Hz.

Test this yourself by feeding the DFT a pure sine wave at a particular frequency, and see if it matches your calculation.

I'm not sure what a good value of threshold would be. I would just do the std::max_element of the magnitudes of the frequencies.

Amplitude would be calculated in the time-domain, it's usually done with some sort of rolling max filter, iirc.