Are my normalizations wrong?

Greetings everyone!

I'm working on a project where i need to compare the output of three metods, the FFT, the Lomb-Scargle (LS) periodogram and the Generalised Lomb-Scargle (GLS). The last one is a new metod for the LS periodogram which compansate for measurementes errors and other stuff.

Now, i have the times and values of the signal from 40 subjets and i have to do a frequency analysis. The first metod is the FFT, i made a for loop which calculate the transform for each set of data assuming a uniform sampling (which is false, that's why a will use the LS),calculate the module squared and normalising by the signal length and frequency to obtain the PSD. The last step is to normalise everything for the DC component. Finally it plots every results on the same axes. So far so good.

t=tDX;

Accomodamento=b1

dt=diff(t);% Delta t

dt(dt<=0)=NaN;

t(t<=0)=NaN;

Accomodamento(t<=0 | Accomodamento==0)=NaN;

L1=sum(~isnan(t),1);

dt=mean(dt, 'omitnan');

fc=1./dt;

%fn=fc/2; %freq di Nyquist

L=2.^nextpow2(L1);

% DFT

figure(1)

figure(2)

for i=1:width(Accomodamento)

Accomodamento1=Accomodamento(:,i);

Accomodamento1=Accomodamento1(~isnan(Accomodamento1));

L2=L(i);

y = fft(Accomodamento1,L2);

f = fc(i)*(0:(L2/2))/L2; % frequency range

Amp1 = abs(y)/(L1(i)); %Amplitude

%Power1 = (1/((fc(i))*(L1(i)))) * abs(y).^2; % power density

Power1 = (1/(L1(i))) * abs(y).^2;

Power1=Power1./(2*var(Accomodamento1)); %tries to normalise by the variance like the plomb and GLS

Power1=Power1/Power1(1,:);

%single sided DFT plot

Amp=Amp1(1:L2/2+1,:);

Amp(2:end-1,:) = 2*Amp(2:end-1,:);

Power=Power1(1:L2/2+1);

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

figure(1)

hold on

stem(f,Amp);

xlabel('Frequency (Hz)','FontWeight', 'bold')

ylabel('Amplitude(D) ','FontWeight', 'bold')

set(gca,'FontWeight','bold')

title("Right eye with FFT " + range0 +"D-" + range1 + "D")

xlim([0 5]);

hold off

figure(2)

hold on

stem(f(2:end),Power(2:end))

xlabel('Frequency (Hz)','FontWeight', 'bold')

%ylabel('PSD(D^{2}/Hz)','FontWeight', 'bold')

ylabel('Normalized PSD','FontWeight', 'bold')

set(gca,'FontWeight','bold')

title("Right eye with FFT " + range0 +"D-" + range1 + "D")

xlim([0 5]);

hold off

end

figure(1)

TFplot1=gca;

figure(2)

TFplot2=gca;

Now i want to use the plom function, i tried the 'psd' option and the 'normalized' option, and the results are way too different from the FFT (normalized of not), concerning the values on the Y axes. i know it may be because the LS accounts for the uneven sampling interval, but i expected a similar results on the y axes. Her's the question, since the plomb function does not calculate the DC value (which i made manually as the mean of the signal squared for the PSD), i can't check if dividing by it can fix the issue.

t=tDX;

Accomodamento=b1

t(t<=0)=NaN;

Accomodamento(t<=0 | Accomodamento==0)=NaN;

Accomodamento = Accomodamento(:,~all(isnan(Accomodamento))); %remove NaN col

t = t(:,~all(isnan(t)));

figure(20)

Lombplot2=gca;

hold on

xlim([0 3])

xlabel('Frequency (Hz)','FontWeight', 'bold')

ylabel('Normalized PSD','FontWeight', 'bold')

%ylabel('PSD(D^{2}/Hz)','FontWeight', 'bold')

title("Right eye PSD with Lomb-Scargle" + range0 +"D-" + range1 + "D")

figure(21)

Lombplot1=gca;

hold on

xlim([0 3])

xlabel('Frequency (Hz)','FontWeight', 'bold')

ylabel('Amplitude','FontWeight', 'bold')

title("Right eye Amplitude with Lomb-Scargle" + range0 +"D-" + range1 + "D")

for i=1:width(Accomodamento)

AccomodamentoLomb=Accomodamento(:,i);

AccomodamentoLomb=AccomodamentoLomb(~isnan(AccomodamentoLomb));

tLomb=t(:,i);

tLomb=tLomb(~isnan(tLomb));

[yLomb, fLomb]=plomb(AccomodamentoLomb,tLomb,'normalized');

% AmpLomb=sqrt(yLomb.*fLomb);

AmpLomb=sqrt(yLomb.*var(AccomodamentoLomb));

figure(20)

stem(fLomb,yLomb)

figure(21)

stem(fLomb,AmpLomb)

clear AccomodamentoLomb tLomb yLomb fLomb

end

hold off

clear t Accomodamento

Lastly i tried using the code for the GLS, made by the same researcher the introduced it (from THIS paper). Aside from a problem whth the time vector, esaly fixed, it works perfectly, but again the y values are too different from the FFT ones. the code i used is the same ad the plomb one but it calls the code below instead of the plomb function

function [f,pow,AMP,fap] = gls_MZ(t, y, f, e, fap_vals)

%

% Input: t - time vecotor of the data

% y - the data

% f - frequencies vector to calculate GLS (optional)

% e - errors of the data (optional)

% fap_vals - GLS power values to calculate the fap (default = [0.1 0.01 0.001])

%

% Output: f - frequencies.

% pow - normalized GLS power for each f.

% fap - selected false-alarm probability GLS values.

%

% Created: 2015, Shay Zucker

% Modified: 2018, Lev Tal-Or & Mathias Zechmeister

ofac = 2;

hifac = 1;

fstep = 1 / tbase / ofac ; % frequency sampling depends on the time span, default for start frequency

%fbeg = fstep;

fbeg = 0;

fnyq = 0.5 / tbase * N; % Nyquist frequency

fend = fnyq * hifac;

f = fbeg:fstep:fend;

% Allocate output array.

%

pow = 0 * f;

AMP = 0*f;

% Convert frequencies to angular frequencies.

%

w = 2 * pi * f;

Y = mean(y);

YY_hat = mean(y.^2);

YY = YY_hat - Y^2;

%YY = var(y);

for k = 1:length(f)

x = w(k) * t; % phases

cosx = cos(x);

sinx = sin(x);

cos2x = cos(2*x);

sin2x = sin(2*x);

omegatau = atan2(mean(sin2x)-2*mean(cosx)*mean(sinx),...

mean(cos2x)-(mean(cosx)^2-mean(sinx)^2))/2;

cosx_tau = cos(x-omegatau);

sinx_tau = sin(x-omegatau);

C = mean(cosx_tau);

S = mean(sinx_tau);

YC_hat = mean(y.*cosx_tau);

YS_hat = mean(y.*sinx_tau);

YC = YC_hat - Y*C;

YS = YS_hat - Y*S;

CC_hat = mean(cosx_tau.^2);

SS_hat = mean(sinx_tau.^2);

CC = CC_hat - C^2;

SS = SS_hat - S^2;

pow(k) = (1/YY * (YC^2/CC+YS^2/SS));

%my extra lines below

powDC= (Y^2) / YY; %all normalised by variance

%pow(k)= pow(k)/powDC;

AMP(k)=sqrt(pow(k)*YY*(Y^2));

end

I'd like to know which are the relations between the FFT normalization and the plomb and GLS ones, cause for this study the power and amplitude values are important. I tried to normalize every output whith the same values but it's even worse. I'm clearly missing something but i don't know what i can try now

I'll attach time and signal data for example.

Thank you in advance!

[Jan]

It increases the chance to get a useful answer, if a reader has the chance to understand, what you are asking for, in less than a minute.

You mentioned some details, which can be understood only by persons, who work exactly in your field. As far as I understand, your problem has no relation to Matlab, but it a physical question?

[Andrea Carobbi]

Yes and no actually. I mean, it's a question born from a physical problem but my main concern is if i'm missing something in the code and the way the functions work. For example, if i remove all the normalizations in the FFT and in plomb (to talk only about MATLAB functions) and calculate the PSD as reported in the help documentation, the results on the Y axes are very different (which was expected, but not as much as they are). I'm cleary wrong on how this functions work and i'm missing something, and the MATLAB documentation doesen't help the way i hoped.

I'm really sorry the question was too dense, i tried writing any information that might be usefull without going too deep.

Sometimes i found questions similarly written that had exactly the answer i was looking for, so i hoped this could be the case.

[dpb]

Well, there are about as many ways to normalize an FFT as there are users -- that's an exaggeration of course, but it all depends on what your application is and what quantity you're trying to measure.

The example PSD in the MATLAT documentation normalizes the one-sided PSD such that the bin peaks match the amplitude in the time domain -- they use a 0.7 and 1.0 peak sine wave and for the pure signal the result returned is exactly (well, within double precision accuracy) 0.7 and 1.0.

How that compares to what is done in those other methods you'll have to dig into their particular features to know, but that isn't typical for power measurements; it's a valid normalization, but it isn't directly related to power, but to the input waveform amplitide.

Also note the effects of noise and sampling frequency and signal sample length -- unless the input frequencies are identially the same as one of the FFT bins (as they were chosen to be in the example), the energy/amplitude of the FFT will be spread out into the adjacent bins and so the total amplitude of a given peak is the integral of the peak, not just the single maximum point/bin value. And, as the example illustrates, adding noise changes the comparison such that the values are not identical to the noise-free amplitudes. One can improve this with averaging, but it's a difference.

I've written innumerable Answers here; many of which illustrate the above points--although they're easy-enough to see simply by changing some parameters in the example -- change the sine input frequency by 1 Hz, say, and rerun and see the difference.

[Andrea Carobbi]

Sorry for the late reply and thanks for the answer.

I know the normalization doesn't have a "standard" and that's what is bothering me. Honestly i tried to normalize because the plots was extremly different from what i obtained whith other metods, so i though "maybe if i normalize everthing whith the same factors at least i can compare them". Nope, wrong scale again, so i though i was wrong in using the function or the function was doing something "strange" i coundn't understand. More so i normalized by the DC value in FFT but plomb can't compute it (it gets a 0/0 at 0Hz), so i'm even more confused...

But anyway, if the code isn't wrong and i did everything correctly, like Jan said this is not a MATLAB question anymore, so thank you very much for answering me. I submitted the results i got and waiting for the feedback, there's no much more i can do.

[dpb]

"...so i normalized by the DC value in FFT but plomb can't compute it (it gets a 0/0 at 0Hz), ..."

Subtract the mean from the time series first and you'll get zero mean from the FFT -- but doing so won't change the magnitude of the other frequency components at all; you can just as easily simply zero-out the DC bin in the computed FFT.

[Andrea Carobbi]

Uh, thank you! I tried doing this on the FFT and leaving averything else as it is and now the plots are much better, they're at least on the same order of magnitude. So this is what i was doing wrong...