From: glowkeeper on
The planning function I'm using is fftwf_plan_r2r_1d, and I use it this
way:

fftwf_plan_r2r_1d( 2048, m_pFourierIn, m_pFourierOut, FFTW_R2HC,
FFTW_BACKWARD );

But nonetheless, I'm still getting the power algorithm wrong it
seems....

>>I suppose I could also do something like taking the square root of the
>>square, but that seems expensive (although granted, I don't know what
>fabs
>>is doing under the hood).
>
>First of all, your Fourier coefficients should be complex, even if the
>input data are real. From FFTW's website, the planning function is:
>
>fftw_plan fftw_plan_dft_r2c_1d(int n, double *in, fftw_complex *out,
> unsigned flags);
>
>where *in is a pointer to your real input array, and *out is a pointer
to
>your complex output array.
>
>Next, fabs() performs an absolute value of a real data point. For
power,
>you actually want the absolute value squared. I.e., given a complex
>coefficient of (a+jb), the power is (a+jb)(a-jb), or abs(a+jb)^2. Since
>you're using FFTW, the output array is defined as type fftw_complex (or
>fftwf_complex for single-precision), unless you've included the
complex.h
>header from the C99 math library. If so, then it will be of type
>double/float complex at the output. I've never done it this way but
maybe
>Steven Johnson can chime in on how it works better (I prefer their
>convention since the complex math doesn't work as well on the processors
>I've used). Either way, the appropriate function at that point would be
>cabs(x)*cabs(x) (well, y = cabs(x), then y2 = y*y, if you want to
out-think
>the compiler), which is the complex magnitude squared.
>
>Mark
>
From: glowkeeper on
Actually, I should say that I'm using

fftwf_plan_r2r_1d( 2048, m_pFourierIn, m_pFourierOut, FFTW_R2HC,
FFTW_BACKWARD );

and then when I come to do the power computations I work on the first 1024
entries in the returned transform (ignoring the imaginary parts). Is this
incorrect?

>The planning function I'm using is fftwf_plan_r2r_1d, and I use it this
>way:
>
>fftwf_plan_r2r_1d( 2048, m_pFourierIn, m_pFourierOut, FFTW_R2HC,
>FFTW_BACKWARD );
>
>But nonetheless, I'm still getting the power algorithm wrong it
>seems....
>
>>>I suppose I could also do something like taking the square root of the
>>>square, but that seems expensive (although granted, I don't know what
>>fabs
>>>is doing under the hood).
>>
>>First of all, your Fourier coefficients should be complex, even if the
>>input data are real. From FFTW's website, the planning function is:
>>
>>fftw_plan fftw_plan_dft_r2c_1d(int n, double *in, fftw_complex *out,
>> unsigned flags);
>>
>>where *in is a pointer to your real input array, and *out is a pointer
>to
>>your complex output array.
>>
>>Next, fabs() performs an absolute value of a real data point. For
>power,
>>you actually want the absolute value squared. I.e., given a complex
>>coefficient of (a+jb), the power is (a+jb)(a-jb), or abs(a+jb)^2.
Since
>>you're using FFTW, the output array is defined as type fftw_complex (or
>>fftwf_complex for single-precision), unless you've included the
>complex.h
>>header from the C99 math library. If so, then it will be of type
>>double/float complex at the output. I've never done it this way but
>maybe
>>Steven Johnson can chime in on how it works better (I prefer their
>>convention since the complex math doesn't work as well on the
processors
>>I've used). Either way, the appropriate function at that point would
be
>>cabs(x)*cabs(x) (well, y = cabs(x), then y2 = y*y, if you want to
>out-think
>>the compiler), which is the complex magnitude squared.
>>
>>Mark
>>
>
From: markt on
>Actually, I should say that I'm using
>
>fftwf_plan_r2r_1d( 2048, m_pFourierIn, m_pFourierOut, FFTW_R2HC,
>FFTW_BACKWARD );

Ah, gotcha. OK, the way I read the documentation on this is that the
complex data still exists, but in the second half of the transform output,
which you would need to incorporate into your absolute value function, i.e.
you still get a complex output, but the data are ordered all of the inphase
first, then quadrature last. I'll defer to SJ for a better explanation
since I've never used the R2HC option.

Mark
From: Ron N. on
On Apr 24, 8:38 am, "glowkeeper" <glowkee...(a)yahoo.co.uk> wrote:
> Actually, I should say that I'm using
>
> fftwf_plan_r2r_1d( 2048, m_pFourierIn, m_pFourierOut, FFTW_R2HC,
> FFTW_BACKWARD );
>
> and then when I come to do the power computations I work on the first 1024
> entries in the returned transform (ignoring the imaginary parts). Is this
> incorrect?

Yes. If you ignore the "imaginary" parts, then a phase
shift of 90 degrees could cause the fundamental or any
harmonic to disappear from your composite "spectrum".
If your fft frames are not phase locked to the signal
generator, then this could easily cause the random type
of results you may be seeing.



From: Ron N. on
On Apr 23, 11:33 pm, "glowkeeper" <glowkee...(a)yahoo.co.uk> wrote:
> >If you downsample your fft, shouldn't that change the bin
> >width? Or are you throwing away spectral energy? That
> >loss of information could certainly effect the results.
>
> Ah, now this was something I wasn't sure of to be honest. Maybe? I output
> my product to an array the same size as the fft, and therefore I thought
> the bin width was the same.
>
> But I think you're suggestingt that my 2x bin width is ~5Hz/2, 3x ~5Hz/3,
> 4x ~5Hz/4, ...

Actually, the opposite. When you downsample a spectrum,
each bin should end up representing a wider frequency band,
e.g. 5 Hz * 2, etc., or you will be throwing away energy.


IMHO. YMMV.