Calculating magnitude of FFT result

Hi, I hope this is the correct forum for such a question.

I am writing some DSP code that runs on a 32 bit CPU for the purpose of
analysing vibration in a machine. I am acquiring samples from an
accelerometer at a sampling rate of 4KHz using a 16bit ADC. I scale the
values so that the full scale value (32767) from the ADC equals 1.0 in Q20
format. I then apply the Hamming window and calculate a 1024 point FFT.
Using the results from the FFT I examine the frequencies of interest using
the quadratic method to interpolate the frequency and magnitude allowing
for the coherent gain of the window. While I get an accuracy of 1Hz for the
frequency the magnitude is less than expected due to the Hamming window.
While I have read that the equivalent window function can be applied in the
frequency domain instead of the time domain I haven't quite been able to
work out the math. One other solution I have come up with is to calibrate
the system at power up by simulating a sine wave of a specific amplitude at
each frequency (half of the frequencies due to symmetry) and calculating
the magnitude correction multipliers and store them in a lookup table. Is
this a reasonable approach or am I overlooking something?

Thanks

	 

_____________________________		
Posted through www.DSPRelated.com

On Thu, 02 Oct 2014 09:30:40 -0500, "Robd" <101885@dsprelated> wrote:

>Hi, I hope this is the correct forum for such a question.
>
>I am writing some DSP code that runs on a 32 bit CPU for the purpose of
>analysing vibration in a machine. I am acquiring samples from an
>accelerometer at a sampling rate of 4KHz using a 16bit ADC. I scale the
>values so that the full scale value (32767) from the ADC equals 1.0 in Q20
>format. I then apply the Hamming window and calculate a 1024 point FFT.
>Using the results from the FFT I examine the frequencies of interest using
>the quadratic method to interpolate the frequency and magnitude allowing
>for the coherent gain of the window. While I get an accuracy of 1Hz for the
>frequency the magnitude is less than expected due to the Hamming window.
>While I have read that the equivalent window function can be applied in the
>frequency domain instead of the time domain I haven't quite been able to
>work out the math. One other solution I have come up with is to calibrate
>the system at power up by simulating a sine wave of a specific amplitude at
>each frequency (half of the frequencies due to symmetry) and calculating
>the magnitude correction multipliers and store them in a lookup table. Is
>this a reasonable approach or am I overlooking something?
>
>Thanks

Sounds reasonable to me.  I've done similar things that worked well
for their intended application.   I'm always a big fan of calibration
since it can potentially compensate for known and unknown distortions.


Eric Jacobsen
Anchor Hill Communications
http://www.anchorhill.com

Multiplication in time domain <=> circular convolution in frequency domain
In Matlab:
xr = complex(randn(1024,1), randn(1024,1));
win = hamming(1024,'periodic');
X1 = fft(win.*xr);
X2 = fft(xr);
X2w = filter([-0.23 0.54 -0.23],1, [X2(end-1:end);X2;X2(1:2)]); % Hamming win
plot(([X1 X2w(3+(1:1024))]),'.-');
---
rgds

On Thu, 02 Oct 2014 09:30:40 -0500, Robd wrote:

> Hi, I hope this is the correct forum for such a question.
> 
> I am writing some DSP code that runs on a 32 bit CPU for the purpose of
> analysing vibration in a machine. I am acquiring samples from an
> accelerometer at a sampling rate of 4KHz using a 16bit ADC. I scale the
> values so that the full scale value (32767) from the ADC equals 1.0 in
> Q20 format. I then apply the Hamming window and calculate a 1024 point
> FFT. Using the results from the FFT I examine the frequencies of
> interest using the quadratic method to interpolate the frequency and
> magnitude allowing for the coherent gain of the window. While I get an
> accuracy of 1Hz for the frequency the magnitude is less than expected
> due to the Hamming window. While I have read that the equivalent window
> function can be applied in the frequency domain instead of the time
> domain I haven't quite been able to work out the math. One other
> solution I have come up with is to calibrate the system at power up by
> simulating a sine wave of a specific amplitude at each frequency (half
> of the frequencies due to symmetry) and calculating the magnitude
> correction multipliers and store them in a lookup table. Is this a
> reasonable approach or am I overlooking something?

If by "simulated" you just mean that you're doing it all numerically, then 
as long as you don't change the window, you could just do it at build time 
or design time.

For a simple enough window you could also derive a symbolic expression for 
the attenuation and use that, although from a practical perspective that 
may not be the best use of your time.

-- 

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com

On Thu, 02 Oct 2014 09:30:40 -0500, "Robd"
<101885@dsprelated> wrote:

>Hi, I hope this is the correct forum for such a question.
>
>I am writing some DSP code that runs on a 32 bit CPU for the purpose of
>analysing vibration in a machine. I am acquiring samples from an
>accelerometer at a sampling rate of 4KHz using a 16bit ADC. I scale the
>values so that the full scale value (32767) from the ADC equals 1.0 in Q20
>format. I then apply the Hamming window and calculate a 1024 point FFT.
>Using the results from the FFT I examine the frequencies of interest using
>the quadratic method to interpolate the frequency and magnitude allowing
>for the coherent gain of the window. While I get an accuracy of 1Hz for the
>frequency the magnitude is less than expected due to the Hamming window.
>While I have read that the equivalent window function can be applied in the
>frequency domain instead of the time domain I haven't quite been able to
>work out the math. One other solution I have come up with is to calibrate
>the system at power up by simulating a sine wave of a specific amplitude at
>each frequency (half of the frequencies due to symmetry) and calculating
>the magnitude correction multipliers and store them in a lookup table. Is
>this a reasonable approach or am I overlooking something?
>
>Thanks
>

The time domain is the simplest place to apply the window is
in the time domain... it's just a sample-by-sample multiply
there.

The window gain change is related to the area under the
window.  You can probably find this tabulated somewhere for
standard windows like Hann, Hamming, etc.  If you create and
store the N window values in a table ahead of time, it's
easy to get the integral then. You shouldn't need to do your
simulation/calibration for this, since the window function
is known.

Instead of the quadratic method to determine peak frequency
and amplitude, I prefer a table-based interpolation method.
Essentially, I did what you describe for your calibration,
stepping the frequency in fractional-line increments and
noting the ratio of the peak amplitude and the highest
adjacent line.  You only need to do this for one span of
frequencies between adjacent spectrum lines, but you need to
do this mid-spectrum... if you are too close to DC you get
interference from negative frequencies.  (And that means the
resulting interpolation doesn't work well in the lowest few
spectral lines.)

Best regards,









Bob Masta
 
              DAQARTA  v7.60
   Data AcQuisition And Real-Time Analysis
              www.daqarta.com
Scope, Spectrum, Spectrogram, Sound Level Meter
 Frequency Counter, Pitch Track, Pitch-to-MIDI 
   FREE Signal Generator, DaqMusiq generator    
          Science with your sound card!

On Thursday, October 2, 2014 10:30:40 AM UTC-4, Robd wrote:
> Hi, I hope this is the correct forum for such a question.
> 
> 
> 
> I am writing some DSP code that runs on a 32 bit CPU for the purpose of
> 
> analysing vibration in a machine. I am acquiring samples from an
> 
> accelerometer at a sampling rate of 4KHz using a 16bit ADC. I scale the
> 
> values so that the full scale value (32767) from the ADC equals 1.0 in Q20
> 
> format. I then apply the Hamming window and calculate a 1024 point FFT.
> 
> Using the results from the FFT I examine the frequencies of interest using
> 
> the quadratic method to interpolate the frequency and magnitude allowing
> 
> for the coherent gain of the window. While I get an accuracy of 1Hz for the
> 
> frequency the magnitude is less than expected due to the Hamming window.
> 
> While I have read that the equivalent window function can be applied in the
> 
> frequency domain instead of the time domain I haven't quite been able to
> 
> work out the math. One other solution I have come up with is to calibrate
> 
> the system at power up by simulating a sine wave of a specific amplitude at
> 
> each frequency (half of the frequencies due to symmetry) and calculating
> 
> the magnitude correction multipliers and store them in a lookup table. Is
> 
> this a reasonable approach or am I overlooking something?
> 
> 
> 
> Thanks
> 
> 
> 
> 	 
> 
> 
> 
> _____________________________		
> 
> Posted through www.DSPRelated.com

Robd,

The math gets easy if you use the definition of the window

w(n)=0.54+0.46*cos(2*pi*n/N), for 0<=n<N, N=fft size.

Check your window and see if the second and last point are the same, if not you are probably using 

w(n)=0.54+0.46*cos(2*pi*n/(N-1)), for 0<=n<N

With the first definition, the fft of the window falls into only 3 bins.

If you are using MATLAB to generate the window, use the 'periodic' option in hamming(). See examples below.

Dirk

>> hamming(7)
ans =
    0.0800
    0.3100
    0.7700
    1.0000
    0.7700
    0.3100
    0.0800
>> fft(hanning(7))/7
ans =
   0.5714          
  -0.2031 - 0.0978i
  -0.0087 - 0.0110i
  -0.0006 - 0.0028i
  -0.0006 + 0.0028i
  -0.0087 + 0.0110i
  -0.2031 + 0.0978i
>> hamming(7,'periodic')
ans =
    0.0800
    0.2532
    0.6424
    0.9544
    0.9544
    0.6424
    0.2532
>> fft(hamming(7,'periodic'))/7
ans =
    0.5400
   -0.2300
         0
         0
         0
         0
   -0.2300
>> 

Thank you for all of the replies!

Yes, by "simulating" I meant numerically. At this point I am just trying to
nullify the effect of the window and get a more flat response from the FFT.
This will hopefully allow a more meaningful calibration when applying
frequencies to the AFE and measuring the response of the anti-aliasing
filter etc.

I was using the ...n)/(N-1) version of the window function as this is what
I found in most of the published material I have read. Since changing to
the ...n)/N versions I have found the effects of the window to be much more
confined (about 3 bins on either side of the main peak). This seems to be
working much better.

Yes, I could calculate the correction multipliers once and either embed
them in the code or save them to flash or EEPROM on first run which I will
probably do eventually.

Does anyone here have any experience with vibration monitoring? More
specifically vibration due to bearing failure? I am currently monitoring
the magnitudes at frequencies determined by the shaft RPM and the bearing
multipliers. In total 6 multipliers plus the shaft speed. The shaft speed
varies to a maximum of 4000 RPM and the highest bearing multiplier is 8.33.
So my max frequency should be around 555Hz. I am sampling at 4KHz and plan
on using a 1024 point FFT. This should allow me to detect frequencies as
low as 4Hz and as high as 1996Hz. Any advise?

Thanks,
Rob
	 

_____________________________		
Posted through www.DSPRelated.com

On Mon, 06 Oct 2014 10:29:26 -0500, "Robd" <101885@dsprelated> wrote:

>Thank you for all of the replies!
>
>Yes, by "simulating" I meant numerically. At this point I am just trying to
>nullify the effect of the window and get a more flat response from the FFT.
>This will hopefully allow a more meaningful calibration when applying
>frequencies to the AFE and measuring the response of the anti-aliasing
>filter etc.
>
>I was using the ...n)/(N-1) version of the window function as this is what
>I found in most of the published material I have read. Since changing to
>the ...n)/N versions I have found the effects of the window to be much more
>confined (about 3 bins on either side of the main peak). This seems to be
>working much better.
>
>Yes, I could calculate the correction multipliers once and either embed
>them in the code or save them to flash or EEPROM on first run which I will
>probably do eventually.
>
>Does anyone here have any experience with vibration monitoring? More
>specifically vibration due to bearing failure? I am currently monitoring
>the magnitudes at frequencies determined by the shaft RPM and the bearing
>multipliers. In total 6 multipliers plus the shaft speed. The shaft speed
>varies to a maximum of 4000 RPM and the highest bearing multiplier is 8.33.
>So my max frequency should be around 555Hz. I am sampling at 4KHz and plan
>on using a 1024 point FFT. This should allow me to detect frequencies as
>low as 4Hz and as high as 1996Hz. Any advise?

Sample more slowly.   If you don't have tight latency constraints, you
can get better resolution from the 1024-pt FFT by reducing the
sampling rate at its input.   Since your max frequency of interest is
only 555 Hz, dropping the sample rate into the FFT by a factor of two
or three will still provide plenty of oversampling and decrease the
bin width by a factor of two or three.   Naturally you'll have to
filter the input signal a bit to prevent aliasing in the sample rate
reduction, but that is doable as a decimating pre-process step before
the FFT.

You can also estimate frequency to within a small fraction of a
bin-width using various interpolation techniques if that is useful to
you.



>Thanks,
>Rob
>	 
>
>_____________________________		
>Posted through www.DSPRelated.com

Eric Jacobsen
Anchor Hill Communications
http://www.anchorhill.com

On Mon, 06 Oct 2014 16:59:18 GMT, eric.jacobsen@ieee.org
(Eric Jacobsen) wrote:

>On Mon, 06 Oct 2014 10:29:26 -0500, "Robd" <101885@dsprelated> wrote:
>
>>Thank you for all of the replies!
>>
>>Yes, by "simulating" I meant numerically. At this point I am just trying to
>>nullify the effect of the window and get a more flat response from the FFT.
>>This will hopefully allow a more meaningful calibration when applying
>>frequencies to the AFE and measuring the response of the anti-aliasing
>>filter etc.
>>
>>I was using the ...n)/(N-1) version of the window function as this is what
>>I found in most of the published material I have read. Since changing to
>>the ...n)/N versions I have found the effects of the window to be much more
>>confined (about 3 bins on either side of the main peak). This seems to be
>>working much better.
>>
>>Yes, I could calculate the correction multipliers once and either embed
>>them in the code or save them to flash or EEPROM on first run which I will
>>probably do eventually.
>>
>>Does anyone here have any experience with vibration monitoring? More
>>specifically vibration due to bearing failure? I am currently monitoring
>>the magnitudes at frequencies determined by the shaft RPM and the bearing
>>multipliers. In total 6 multipliers plus the shaft speed. The shaft speed
>>varies to a maximum of 4000 RPM and the highest bearing multiplier is 8.33.
>>So my max frequency should be around 555Hz. I am sampling at 4KHz and plan
>>on using a 1024 point FFT. This should allow me to detect frequencies as
>>low as 4Hz and as high as 1996Hz. Any advise?
>
>Sample more slowly.   If you don't have tight latency constraints, you
>can get better resolution from the 1024-pt FFT by reducing the
>sampling rate at its input.   Since your max frequency of interest is
>only 555 Hz, dropping the sample rate into the FFT by a factor of two
>or three will still provide plenty of oversampling and decrease the
>bin width by a factor of two or three.   Naturally you'll have to
>filter the input signal a bit to prevent aliasing in the sample rate
>reduction, but that is doable as a decimating pre-process step before
>the FFT.
>
>You can also estimate frequency to within a small fraction of a
>bin-width using various interpolation techniques if that is useful to
>you.

The OP may want to consider using an ordinary Windows sound
card.   Sound cards are AC coupled, but it's easy to find
one that goes to 4 Hz.  See my tests of several USB devices
at <http://www.daqarta.com/dw_gguu.htm>.  One of them (7.1
channel CM6206)  goes below 1 Hz, and it's only US$20.

The default sample rate for Windows sound is 48000 Hz, but
if you set it lower Windows will take care of the sample
rate conversion and anti-aliasing automagically.  I'd
recommend an integer submultiple of 48000, such as 2400.

Shameless plug:  My Daqarta software will probably do
everything you need, including peak frequency estimation to
about 1/50 bin-width.  You get immediate real-time waveform,
spectrum, or spectrogram displays with lots of triggering
options, cursor readouts, and all the details that can
otherwise suck up a lot of coding time.

Even if you eventually decide you need something it can't do
(either as-installed, or with the built-in macro system),
you may still find that it's a good way to ease into this
project.   You can try it free, with all features, for 30
sessions / 30 days, and I will be glad to renew that upon
request.  The Pro version is only US$99 if you decide not to
roll your own from scratch.

Best regards,


Bob Masta
 
              DAQARTA  v7.60
   Data AcQuisition And Real-Time Analysis
              www.daqarta.com
Scope, Spectrum, Spectrogram, Sound Level Meter
 Frequency Counter, Pitch Track, Pitch-to-MIDI 
   FREE Signal Generator, DaqMusiq generator    
          Science with your sound card!

On Monday, October 6, 2014 11:29:26 AM UTC-4, Robd wrote:
> Thank you for all of the replies!
> 
> 
> 
> Yes, by "simulating" I meant numerically. At this point I am just trying to
> 
> nullify the effect of the window and get a more flat response from the FFT.
> 
> This will hopefully allow a more meaningful calibration when applying
> 
> frequencies to the AFE and measuring the response of the anti-aliasing
> 
> filter etc.
> 
> 
> 
> I was using the ...n)/(N-1) version of the window function as this is what
> 
> I found in most of the published material I have read. Since changing to
> 
> the ...n)/N versions I have found the effects of the window to be much more
> 
> confined (about 3 bins on either side of the main peak). This seems to be
> 
> working much better.
> 
> 
> 
> Yes, I could calculate the correction multipliers once and either embed
> 
> them in the code or save them to flash or EEPROM on first run which I will
> 
> probably do eventually.
> 
> 
> 
> Does anyone here have any experience with vibration monitoring? More
> 
> specifically vibration due to bearing failure? I am currently monitoring
> 
> the magnitudes at frequencies determined by the shaft RPM and the bearing
> 
> multipliers. In total 6 multipliers plus the shaft speed. The shaft speed
> 
> varies to a maximum of 4000 RPM and the highest bearing multiplier is 8.33.
> 
> So my max frequency should be around 555Hz. I am sampling at 4KHz and plan
> 
> on using a 1024 point FFT. This should allow me to detect frequencies as
> 
> low as 4Hz and as high as 1996Hz. Any advise?
> 
> 
> 
> Thanks,
> 
> Rob
> 
> 	 
> 
> 
> 
> _____________________________		
> 
> Posted through www.DSPRelated.com

sounds like you should use flat top windows. no compensation needed.

Clay