How could it be that by FFT is not symmetric?

I perform fft as follows

dt = 1/Fs;
% dt = 0.216;
N = numel(rho{icell,iday});
df = 1/(dt*N);
FTrho = fft(rho{icell,iday},N);
powT = sum(abs(rho{icell,iday}).^2)*dt
powF = sum(abs(FTrho*dt).^2)*df

powT =

     4.5455e+06

powF =

     4.5455e+06

the last two lines is just to check that the energy is preserved. when I plot:

figure;plot(abs(FTrho),'-r');xlim([-400,1500])

I get this:

when I do as follows:

FTrho(2:end) = FTrho(2:end)*2;

my energy in time and frequency domain are not equal:

powT =

     4.5455e+06

powF =

     6.7784e+06

I want to claculate the energy of certains fft peaks, therefore I want to be sure my fft yaxis is correctly scalled.

thank you very much for your,

Tamara

[Matt J]

Since you already had correct scaling (the energies were equal) to begin with, why would you disturb that?

[topolowa]

Hi Matt,

I want to calculate energy of certain peaks in fft, not all the spectrum. The energy of given frequency peak has an exact physical meaning, therefore I want to set my yaxis scale right. I understand your comment, normally I would care about a spectrum shape, rather than its absolute values. However not this time. Normally I would take halve of FFT, but which halve should I take if they are not equal?

thanks

Tamara

Answer 1

[dpb]

FTrho(2:end-1) = 2*FTrho(2:end-1);

You're doubling the DC component, too.

Answer 2

[Matt J]

The quantity

FTrho = fft(rho{icell,iday},N)*dt;

is a correctly-scaled approximation to the continuous Fourier transform. The energy of particular regions of the spectrum would be computed by integrating over that region.

[topolowa]

Hi Guys, thanks for you time.

would you then agree with the way I calculate the energy (see code below)?

for icell = 1:ncell
    for iday = 1:nday
        tmpfft = fft(rho{icell,iday},N);
        FFTrho{icell,iday} = tmpfft(1:floor(N/2));
        FFTrhoAbs{icell,iday} = abs(Ts*FFTrho{icell,iday}(1:end));
        FFTrhoAbs{icell,iday}(2:end) = FFTrhoAbs{icell,iday}(2:end)*2;
        f = (0:N-1)*Fs/N;
        xfft = f(1:numel(FFTrho{icell,iday}));
        % find harmonics
        [pks{icell,iday}, locs{icell,iday},pksw{icell,iday}, pksprom{icell,iday}] = ...
        findpeaks((FFTrhoAbs{icell,iday}),xfft,'MinPeakProminence',promin,'MinPeakHeight',minH);

          absFFT2 = FFTrhoAbs{icell,iday}.^2;
          df = Fs/N;
          % calculate energy of given peak/harmonic
          for ipeak = 1:npeak
              energy{icell,iday}(ipeak,1) = trapz(xfft(indst:inden),df*absFFT2(indst:inden));
          end
      end
  end

where indst, inden are start and end indices for a given peak.

cheers

Tamara

[Matt J]

Seems right to me.

[dpb]

It is, however, 2-sided so half the energy is in the negative frequency range and half positive. If want peaks on 0-Fmax to be full magnitude then do need the "times-two" factor.

There's example in the FFT doc that illustrates getting back the identical peak magnitude for noiseless signal that precisely matches binning; noise and frequency smearing affect that in real world, of course.

[topolowa]

Hi, I use one sided FFT:

tmpfft(1:floor(N/2))

and only after, multiply by 2,

or I am wrong...maybe I should use fftshift first. thanks

Tamara

[dpb]

See earlier Answer; the FFT is two-sided but there's only 1 DC component; not two so doubling 1:L/2 doubles it twice as it is already full magnitude which is where the initial asymmetry came from. If there is a significant DC component in the input signal, then that doubling is enough to be apparent as in your case; if it were tiny or took the DC component out before applying FFT the effect would likely not be noticed.

From the signal in the doc example, look at P2=abs(S/L) --

>> [P2(1:10);P2(end:-1:end-9)]
ans =
   1.0e-15 *
0.0600  0.0372  0.3311   0.3308   0.1617   0.2704   0.0999   0.1121  0.0886    0.2677
0.0372  0.3311  0.3308   0.1617   0.2704   0.0999   0.1121   0.0886   0.2677   0.1087
>>

NB: the symmetry after the first element (DC) to the end; there's no DC component at the end of the returned FFT vector.

Similarly for Fmax the symmetry in the complex FFT is about the L/2+1 point, not L/2 and so:

>> P2((L/2+1)-5:(L/2+1)+5)
ans =
   1.0e-15 *
0.3544 0.3263  0.1500  0.6822  0.2950  0.4871  0.2950  0.6822  0.1500  0.3263  0.3544
>> 

NB: the Fmax element at L/2+1 is unique; not two of 'em and note its magnitude compared to each side.