Re: How do YOU handle this disparity in fft?

[Mike Miller]

Hi, great topic, I'll try to comment briefly on your discussion as you
develop your points, answered inline below.

On Mon, Feb 10, 2014 at 19:07:49 -0800, Macy wrote:
> I'm using fft and want to model everything. To do that, everything needs to 
> be consistent else all falls apart.
> 
> My 'difficulty' is making apples compare to apples. Let's assume a system 
> that acquires data at the rate of 100kS/s, where we'll only use smaller 
> chunks of that data, streaming in at ten times a second, 10000 samples each 
> chunk.
> 
> These are always real waveforms so the spectrum will always be 'symmetrical' 
> about the Nyquist rate of the fft, but there are two places in the fft cause 
> me a problem, one at DC and the other at the Nyquist rate [EVEN sample 
> packets].

Right, because these two FFT bins are the midpoints of the circular
symmetry of the DFT. More further on.

> What I mean is if I adjust everything to make the 'tones' rms, and then 
> appear correct in the fft, while only using the lower half the spectrum; the 
> DC value is wrong. And the value at exatly the Nyquist rate, 
> samplepacketlength/2, is wrong:
> 
> Do the exercise for yourself and see:
> 
> Assume 100kS/s ADC with packet lengths of 10000. 
> 
> To make a 1 Vrms 100Hz tone starting at t=0 is:
> Sr=100000; p=10000;
> t=([1:p]-1)/Sr;
> f=100;
> sig=sqrt(2)*cos(2*pi()*f*t);
> %real wave so need only lower half of solution:
> b=sqrt(2)*fft(sig)/p;
> % note the resolution is 10Hz per index, so index of 100 is 10 + 1
> b(11)
> 
> the answer you get is 1, which is indeed the rms value of the 100Hz tone
> you will ALWAYS get that answer no matter what the frequency of the tone is. 
> 
> But now try that with a DC value of 1 and the answer for b(0) becomes 1.414, 
> not right! and if you happen to have any tones at the Nyquist rate of 5kHz, 
> then at 5001, that 1 Vrms 'looks' like 2

You mean 50 kHz and bin 5001 (10 Hz per bin), right?

The best way to look at this is to recall that you are discarding half
of the signal power when you only look at the one-sided FFT. So for
every FFT bin from k=2:5000, only half of the signal power is present.
Your multiplying by sqrt(2) is correcting for the other half of the
power in the bins from k=5002:10000. But this obviously shouldn't apply
to the DC bin because it is complete in itself.

As for bin 5001, its power is double what you expect because it contains
the total energy for that frequency component rather than being split
between two symmetric positive and negative frequency components. It is
its own symmetry point, it contains all the energy from the 50 kHz
sinusoid instead of half, and your normalization results in twice the
power.

>From a practical standpoint, sampling a signal that is exactly Fs/2
gives completely ambiguous results and is best avoided by filtering,
unless you happen to know and be locked exactly to the phase of the
signal.

> So far I have ignored the Nyquist rate because I've just assumed some amount 
> of low pass filtering. And the DC has been no problem, because the DC has 
> been zero, everything is either AC coupled.
> 
> Now how to add Noise Density Function. Which has NO low pass filtering?  
> If I use the NDF*sqrt(Sr/2); all turns out copacetic in appearance.
> note that the average of randn(1,p) is very close to zero and the rms value 
> is 1 
> arbitrarily use density function of 1uV/rtHz 
> 
> sn=sig+1e-6*sqrt(Sr/2)*randn(1,p);
> bn=sqrt(2)*fft(sn)/p;
> % if you do a true rms value on resulting bn noise floor, you get the 
> expected 1uVrms*sqrt(10)
> bnlog=20*log10(abs(bn(1:round(p/2))));
> plot(bnlog)
> % or add the true frequency ordinate, remember the first is DC and is zero 
> freq
> ford=([1:round(p/2)]-1)*Sr/p;
> plot(ford,bnlog);
> 
> The first thing that caught my eye is that even though the total signal plus 
> noise has an average of zero there is 'energy' at DC?  Am I missing something 
> here? 

Only that it's a time-limited sequence, so the average isn't exactly
zero. See below.

> for example, just take a very long random signal, randn(1,1e6), which has an 
> average of zero, no DC then take the energy preserving fft:
> n=sqrt(2)*fft(randn(1,1e6))/1e6;abs(n(1:5))
> the magnitude of n representing the DC value is NOT zero as expected but 
> rather almost equal to the adjacent value(s) ??

Check the mean of your long random signal, it's not really zero and in
fact should be equal to the DC value (scaled by sqrt(2) if you like):

  n = randn (1, 1e6);
  sqrt(2) * mean (n)
  ans =    8.85897299469380e-04
  N = sqrt (2) * fft (n) / 1e6;
  abs (N(1))
  ans =    8.85897299469405e-04

> The next thing is *if* thesignal is DC of 1 Vdc, the fft predicts that it is 
> only .707, not 1 Vrms as expected. 

Not sure what you're doing differently here, but again be sure to scale
sinusoids and not scale DC components by sqrt(2) as appropriate.

> and the last thing is if the signal f is allowed to go to 5kHz it suddenly 
> looks like 2 Vrms, again wrong.

50 kHz again, right?

> However if the tone is anywhere between 10Hz and 4990Hz, the value appears 
> EXACTLY equal to 1 Vrms as expected, but oddly there is NO 'bump' in the 
> middle from the noise through 5kHz.  plot spectrum is flat through the 5kHz 
> and at DC 

Not sure what you mean here, but as I commented above, in the bins from
k=2:5000 what you're doing is exactly correct.

> So far I've been ignoring all this, but have NO  'automatic' way to get 
> around it. 
> I can kind of understand how noise density can fill in a non-existant DC 
> level, the the rest has me wondering.
> 
> [This was a bit long and may have some errors as I presented the question, so 
> hoepfull you can see through those errors and understand my question.]
> 
> For the DC level I should not multiply by 1/sqrt(2), or should I? 
> 
> For an errant tone that has not been LPF'd out, how to handle that 'bump' in 
> the spectrum at p/2+1 ??

Since the amplitude of a signal at exactly Fs/2 can't be trusted anyway,
I'd say you can choose how you want to deal with it: scale it by 1/2 to
counter your earlier scaling, scale it by some arbitrary factor to act
as a digital frequency-domain filter, ignore it.

I hope this helps and answers some of your questions,

-- 
mike