Author Topic: FFT measurements, DC bin and Power  (Read 2306 times)

0 Members and 1 Guest are viewing this topic.

Offline radiolistenerTopic starter

  • Super Contributor
  • ***
  • Posts: 3282
  • Country: ua
FFT measurements, DC bin and Power
« on: March 29, 2023, 05:59:29 pm »
I needs to implement some measurements on ADC sampled data with FFT (such as dBm power in selected bandwidth, etc), and I stuck with one question.

Let's assume that load R=1. So, the power of a sine waveform signal with a peak amplitude A can be calculated as:

Power = A^2 / 2

but at the same time power of DC signal can be calculated as:

Power = A^2

This is because DC signal is constant, so it has constant power and don't needs to be averaged over full cycle with positive and negative part of the waveform.

I'm not sure, may be I'm wrong, but at a glance it seems that a magnitude of DC bin in the FFT result needs to be corrected in order to get proper power value, or at least it needs to be handled in a different way than other FFT bins. I think that such correction is required because FFT DC bin represents constant signal and other FFT bins represents sine carriers.

I asked chat-gpt and he said that I'm correct and DC bin of FFT needs to be corrected for a proper power measurements. But strange thing that I didn't hear about such correction before, so I'm not sure if such correction is relevant for FFT spectrum measurements.

So, I wonder when such correction is needed and when it is not for a spectrum measurements with FFT? What is the best practice?


To be sure, I'm talking about signal which is represented according to AES17 standard meaning, where signal is symmetric around zero value and minimum negative code is not used. So, full scale sine has 0 dBFS level and full scale square wave has +3.01 dBFS level.

And now I'm not sure what is the level of DC signal with maximum positive amplitude? And what is the power of such DC signal?
« Last Edit: March 29, 2023, 06:12:32 pm by radiolistener »
 

Offline gf

  • Super Contributor
  • ***
  • Posts: 1132
  • Country: de
Re: FFT measurements, DC bin and Power
« Reply #1 on: March 29, 2023, 06:09:20 pm »
Consider that there is only one bin for DC and one for Nyquist, but the power of other frequencies is split into positive and negative frequency bins.

EDIT: The sum of the power of the time domain samples is the same as the sum of the power of all (positive and negative!) frequency bins.

Example:

Code: [Select]
>> N=100;
>> signal=rand(1,N);      % a random real signal
>> spectrum=fft(signal);  % corresponding (complex) FFT spectrum

% power sum of time domain samples
>> sum(signal.**2)
ans =  32.211

% power sum of spectrum
>> sum(spectrum.*conj(spectrum))/N
ans =  32.211

% power sum of spectrum (alternative calculation method)
>> sum(abs(spectrum).**2)/N
ans =  32.211
« Last Edit: March 29, 2023, 07:28:38 pm by gf »
 

Online MasterT

  • Frequent Contributor
  • **
  • Posts: 783
  • Country: ca
Re: FFT measurements, DC bin and Power
« Reply #2 on: March 29, 2023, 10:03:40 pm »
I needs to implement some measurements on ADC sampled data with FFT (such as dBm power in selected bandwidth, etc), and I stuck with one question.

Let's assume that load R=1. So, the power of a sine waveform signal with a peak amplitude A can be calculated as:

Power = A^2 / 2

but at the same time power of DC signal can be calculated as:

Power = A^2


You confused A , for AC signal A is peak voltage this is why division by 2 comes. Shoud be RMS, than in both cases DC or AC formula is the same P = A^2
P = Ap-p ^2 / 2 = (A-rms * sqrt(2)) ^2 / 2 = A-rms ^2.
Correction is not required if all magnitudes inside bines are considered as RMS values.
 

Offline gf

  • Super Contributor
  • ***
  • Posts: 1132
  • Country: de
Re: FFT measurements, DC bin and Power
« Reply #3 on: March 29, 2023, 11:01:56 pm »
So, I wonder when such correction is needed and when it is not for a spectrum measurements with FFT? What is the best practice?

If you discard the negative frequency bins and consider only the positive ones, then you need to double the power of each bin (in order to account for the power in the discarded negative frequency counterparts), but don't double the power of the DC and the Nyquist bins (which exist only once and do not have negative frequency counterparts).
 

Offline radiolistenerTopic starter

  • Super Contributor
  • ***
  • Posts: 3282
  • Country: ua
Re: FFT measurements, DC bin and Power
« Reply #4 on: March 30, 2023, 05:39:38 pm »
Here is MATLAB test code, I added 3 test vectors - sine wave, square wave and DC (two of them are commented):
Code: [Select]
Fs = 128;
N = 128;
t = 0:1/Fs:(N-1)/Fs;

signal = sin(2*pi*2*t);                         % sine wave
%signal = cat(2, ones(1, N/2), -ones(1, N/2));  % square wave
%signal = ones(1, N) * 1;                       % DC wave

spectr = fft(signal);

psignal = sum(signal.^2)/N;
pspectr = sum(abs(spectr/N).^2);

fprintf('psignal = %f = %+.2f dBFS\n', psignal, 20*log10(psignal));
fprintf('pspectr = %f = %+.2f dBFS\n', pspectr, 20*log10(pspectr));

the output for sine wave is:
Code: [Select]
psignal = 0.500000 = -6.02 dBFS
pspectr = 0.500000 = -6.02 dBFS

As we see power in time and freq domain are the same which is ok.
But according to AES17 standard, it should be 0 dBFS for sine wave with max amplitude.

How to fix calcualtions to get correct results for AES17 standard?

If assume that 0.5 is a power of a full scale sine wave, we can try to divide power with 0.5 to get level relative to a full scale sine:
Code: [Select]
psignal = (sum(signal.^2)/N) / 0.5;
pspectr = (sum(abs(spectr/N).^2)) / 0.5;

For sine wave it shows 0 dBFS as expected. For square wave and DC it shows the same value +6.02 dBFS. So it looks that there is no need to use correction for DC bin.

But as I remember AES17 standard mention that square wave with max amplitude has +3.01 dBFS level:
Quote
3.3.1
decibels, full scale
dB FS
amplitude expressed as a level in decibels relative to full-scale amplitude (20 times the common logarithm of the amplitude over the full-scale amplitude)
NOTE The rules of the International System of Units (SI) require that a space appear after the standard symbol dB.

3.4
full-scale signal level
FS
signal amplitude relative to the full-scale amplitude expressed in decibels, full scale or percent, full scale
NOTE Because the definition of full scale is based on a sine wave, it will be possible with squarewave test signals to read as much as + 3,01 dB FS. Square-wave signals at this level are not recommended because tilt or overshoot introduced by any filtering operations will cause clipping of the signal.

so I'm not sure where is mistake.

Update: just notice, that AES17 talking about signal level as amplitude...
so I tried to get sqrt to get amplitude from power:
Code: [Select]
asignal = sqrt(sum(signal.^2)/N);
aspectr = sqrt(sum(abs(spectr/N).^2));

fprintf('asignal = %f = %+.2f dBFS\n', asignal, 20*log10(asignal));
fprintf('aspectr = %f = %+.2f dBFS\n', aspectr, 20*log10(aspectr));

it shows the following result:
Code: [Select]
asignal = 0.707107 = -3.01 dBFS
aspectr = 0.707107 = -3.01 dBFS

So, I divided it with 0.707107 (full scale amplitude):
Code: [Select]
asignal = sqrt(sum(signal.^2)/N) * sqrt(2);
aspectr = sqrt(sum(abs(spectr/N).^2))  * sqrt(2);

and it shows expected result:
Code: [Select]
asignal = 1.000000 = +0.00 dBFS
aspectr = 1.000000 = +0.00 dBFS

for both - square wave and DC it shows:
Code: [Select]
asignal = 1.414214 = +3.01 dBFS
aspectr = 1.414214 = +3.01 dBFS

is this correct?
« Last Edit: March 30, 2023, 06:24:34 pm by radiolistener »
 

Offline gf

  • Super Contributor
  • ***
  • Posts: 1132
  • Country: de
Re: FFT measurements, DC bin and Power
« Reply #5 on: March 30, 2023, 07:30:33 pm »
The problem is here: 20*log10(psignal)

For power ratios, dB is 10*log10(psignal)
For voltage ratios , dB is 20*log10(psignal)
 

Offline radiolistenerTopic starter

  • Super Contributor
  • ***
  • Posts: 3282
  • Country: ua
Re: FFT measurements, DC bin and Power
« Reply #6 on: March 31, 2023, 03:17:38 am »
Found that there is required fft result normalization with 1 / sqrt(N). With this normalization rms value calculated in time and frequency domain will be the same:
Code: [Select]
Fs = 128;
N = 128;
t = 0:1/Fs:(N-1)/Fs;

signal = sin(2*pi*2*t);                        % sine wave
%signal = cat(2, ones(1, N/2), -ones(1, N/2));  % square wave
%signal = ones(1, N) * 1;                       % DC wave

spectr = fft(signal) / sqrt(N);

rmsSignal = sqrt(sum(abs(signal).^2)/N);
rmsSpectr = sqrt(sum(abs(spectr).^2)/N);

fprintf('rmsSignal: %f V = %f W\n', rmsSignal, rmsSignal^2);
fprintf('rmsSpectr: %f V = %f W\n', rmsSpectr, rmsSpectr^2);

for sine wave it shows:
Code: [Select]
rmsSignal: 0.707107 V = 0.500000 W
rmsSpectr: 0.707107 V = 0.500000 W

for square wave and DC it shows:
Code: [Select]
rmsSignal: 1.000000 V = 1.000000 W
rmsSpectr: 1.000000 V = 1.000000 W

which looks good, as expected.

So, now in order to get AES17 signal level in dBFS there is need to divide it with rms amplitude of a full scale sine which is 1/sqrt(2), by swapping the numerator and denominator we get multiplier sqrt(2). And then amplitude relation needs to be converted to dB with 20*log10():
Code: [Select]
fprintf('rmsSignal: %+.2f dBFS\n', 20*log10(rmsSignal*sqrt(2)));
fprintf('rmsSpectr: %+.2f dBFS\n', 20*log10(rmsSpectr*sqrt(2)));

for sine wave it shows:
Code: [Select]
rmsSignal: +0.00 dBFS
rmsSpectr: +0.00 dBFS

for square wave and DC it shows:
Code: [Select]
rmsSignal: +3.01 dBFS
rmsSpectr: +3.01 dBFS

By the way, if we use calculations as power it can eliminate redundant sqrt:
Code: [Select]
Fs = 128;
N = 128;
t = 0:1/Fs:(N-1)/Fs;

signal = sin(2*pi*2*t);                        % sine wave
%signal = cat(2, ones(1, N/2), -ones(1, N/2));  % square wave
%signal = ones(1, N) * 1;                       % DC wave

spectr = fft(signal) / sqrt(N);

pSignal = sum(abs(signal).^2)/N;
pSpectr = sum(abs(spectr).^2)/N;

fprintf('rmsSignal: %f V = %f W\n', sqrt(pSignal), pSignal);
fprintf('rmsSpectr: %f V = %f W\n', sqrt(pSpectr), pSpectr);
fprintf('rmsSignal: %+.2f dBFS\n', 10*log10(pSignal*2));
fprintf('rmsSpectr: %+.2f dBFS\n', 10*log10(pSpectr*2));

Code: [Select]
rmsSignal: 0.707107 V = 0.500000 W
rmsSpectr: 0.707107 V = 0.500000 W
rmsSignal: -0.00 dBFS
rmsSpectr: -0.00 dBFS


If I understand correctly, 1/sqrt(N) normalization for FFT is required here to preserve Parseval-Plancherel or Rayleigh's energy theorem.

correct me please if there is some mistake
« Last Edit: March 31, 2023, 05:27:25 am by radiolistener »
 

Offline radiolistenerTopic starter

  • Super Contributor
  • ***
  • Posts: 3282
  • Country: ua
Re: FFT measurements, DC bin and Power
« Reply #7 on: March 31, 2023, 04:21:05 am »
And now, there is need to calculate power of specific bin in the FFT result. I can do it with the following expression:

Code: [Select]
P = 2*abs(spectr(binNumber)).^2/N;

N is FFT size. I'm using multiplier by 2 because there are two parts of spectrum with negative and positive frequency and since I'm use single side spectrum, I need to multiply bin magnitude by 2.

This expression works ok for non DC bin. For example for sine signal it returns the same signal level of specific bin as a total signal level. But for DC signal it returns +6.02 dBFS at DC bin while total rmsSignal = +3.01 dBFS.

So, it seems that there is need to remove multiply by two for DC bin:
Code: [Select]
p0      = abs(spectr(1)).^2/N;
p1      = 2*abs(spectr(2)).^2/N;
p2      = 2*abs(spectr(3)).^2/N;

Is that correct?


By the way, if we use FFT normalization 1/N, it allows to eliminate divide by N for calculations in frequency domain:
Code: [Select]
spectr = fft(signal) / N;

pSignal = sum(abs(signal).^2)/N;
pSpectr = sum(abs(spectr).^2);

fprintf('rmsSignal: %f V = %f W\n', sqrt(pSignal), pSignal);
fprintf('rmsSpectr: %f V = %f W\n', sqrt(pSpectr), pSpectr);
fprintf('rmsSignal: %+.2f dBFS\n', 10*log10(pSignal*2));
fprintf('rmsSpectr: %+.2f dBFS\n', 10*log10(pSpectr*2));

p0      = abs(spectr(1)).^2;
p1      = 2*abs(spectr(2)).^2;
p2      = 2*abs(spectr(3)).^2;
p3      = 2*abs(spectr(4)).^2;
p4      = 2*abs(spectr(5)).^2;
fprintf('p0      = %f V = %+.2f dBFS\n', sqrt(p0),  10*log10(p0*2));
fprintf('p1      = %f V = %+.2f dBFS\n', sqrt(p1),  10*log10(p1*2));
fprintf('p2      = %f V = %+.2f dBFS\n', sqrt(p2),  10*log10(p2*2));
fprintf('p3      = %f V = %+.2f dBFS\n', sqrt(p3),  10*log10(p3*2));
fprintf('p4      = %f V = %+.2f dBFS\n', sqrt(p4),  10*log10(p4*2));

Code: [Select]
rmsSignal: 0.707107 V = 0.500000 W
rmsSpectr: 0.707107 V = 0.500000 W
rmsSignal: +0.00 dBFS
rmsSpectr: +0.00 dBFS
p0      = 0.000000 V = -345.45 dBFS
p1      = 0.707107 V = +0.00 dBFS
p2      = 0.000000 V = -325.27 dBFS
p3      = 0.000000 V = -326.48 dBFS
p4      = 0.000000 V = -331.67 dBFS
« Last Edit: March 31, 2023, 05:58:08 am by radiolistener »
 

Offline gf

  • Super Contributor
  • ***
  • Posts: 1132
  • Country: de
Re: FFT measurements, DC bin and Power
« Reply #8 on: March 31, 2023, 12:28:19 pm »
So, it seems that there is need to remove multiply by two for DC bin:
Code: [Select]
p0      = abs(spectr(1)).^2/N;
p1      = 2*abs(spectr(2)).^2/N;
p2      = 2*abs(spectr(3)).^2/N;

Is that correct?

Yes.

Quote
By the way, if we use FFT normalization 1/N, it allows to eliminate divide by N for calculations in frequency domain:

At least for the Matlab/Octave implementation of FFT. Note that FFT normalization is implementation-dependent, and other FFT libraries may require a different normalization.
 
The following users thanked this post: radiolistener

Offline gf

  • Super Contributor
  • ***
  • Posts: 1132
  • Country: de
Re: FFT measurements, DC bin and Power
« Reply #9 on: March 31, 2023, 09:24:25 pm »
@radiolistener, I also suggest to read this paper which covers even more aspects.
I wanted to post this link earlier, but I couldn't find it right away.

 
The following users thanked this post: rmozel, iMo, radiolistener, mawyatt

Online mawyatt

  • Super Contributor
  • ***
  • Posts: 3194
  • Country: us
Re: FFT measurements, DC bin and Power
« Reply #10 on: April 01, 2023, 06:32:03 pm »
@ gf,

That's a good reference, thanks for posting :-+

Best,
Curiosity killed the cat, also depleted my wallet!
~Wyatt Labs by Mike~
 

Offline radiolistenerTopic starter

  • Super Contributor
  • ***
  • Posts: 3282
  • Country: ua
Re: FFT measurements, DC bin and Power
« Reply #11 on: April 01, 2023, 08:59:11 pm »
@radiolistener, I also suggest to read this paper which covers even more aspects.

Thanks, it contains a lot of info about different window functions, some of them I even didn't hear before :)
 

Offline MegaVolt

  • Frequent Contributor
  • **
  • Posts: 909
  • Country: by
Re: FFT measurements, DC bin and Power
« Reply #12 on: April 05, 2023, 10:51:37 am »
Can someone who is sure of their numbers post some test data? For different window functions?
At the input, for example, 32 numbers; at the output, 16 amplitudes and 16 phases.
 

Offline radiolistenerTopic starter

  • Super Contributor
  • ***
  • Posts: 3282
  • Country: ua
Re: FFT measurements, DC bin and Power
« Reply #13 on: April 05, 2023, 08:53:58 pm »
Can someone who is sure of their numbers post some test data? For different window functions?
At the input, for example, 32 numbers; at the output, 16 amplitudes and 16 phases.

you can do it in MATLAB/Octave/R.

But I notice that MATLAB implementation looks not the best.

For example I wrote my Kaiser window implementation and tried to compare it with MATLAB one. First I found that my function has some minor errors, so I wrote Taylor series approximation for Bessel function which use many terms, and it looks almost the same as MATLAB, but MATLAB has a little difference for a least significant digits.

Here is example of Kaiser window function for length=16 and beta=38:

MATLAB:
Code: [Select]
4.8344100587269438e-16
7.6303044075880228e-09
6.3285566773886172e-06
0.00055999634422098439
0.013171644335785796
0.11723140612283549
0.46884610584259379
0.91996711231932526
0.91996711231932526
0.46884610584259379
0.11723140612283549
0.013171644335785796
0.00055999634422098439
6.3285566773886172e-06
7.6303044075880228e-09
4.8344100587269438e-16

my implementation with Taylor series approximation of Bessel function which use first 100 terms:
Code: [Select]
4.8344100587269438e-16
7.6303044075880245e-09
6.3285566773886138e-06
0.00055999634422098429
0.013171644335785796
0.11723140612283552
0.46884610584259401
0.91996711231932526
0.91996711231932526
0.46884610584259401
0.11723140612283552
0.013171644335785796
0.00055999634422098429
6.3285566773886138e-06
7.6303044075880245e-09
4.8344100587269438e-16

I like Kaiser window, because it allows to see very good dynamic range on FFT, but note that it make bin bandwidth pretty wide. For example, if you use Kaiser window on a sample length N=256 it leads to a more wide bin bandwidth, each bin with max amplitude will have disturbance for at least 12-13 bins on each side, but on the other hand it allows to see details down to -300 dB.

For example, here is a clean sine 32 Hz on FFT with N=256 with applying Kaiser window with beta=38.
And for comparison, here is also screenshot of the same sine on FFT N=256 with applying Blackman-Harris 4 term window:



« Last Edit: April 05, 2023, 09:31:00 pm by radiolistener »
 
The following users thanked this post: MegaVolt


Share me

Digg  Facebook  SlashDot  Delicious  Technorati  Twitter  Google  Yahoo
Smf