How do I calculate area under the curve(power of the spectrum) in frequency axis

clear all
close all
clc
P0=100;
eps0=8.854e-12;
L=10;
n=1.45;
Omega=1.020663297455751e+11;
GAMMAb=3.590391604102621e+08;
c=2.9979e8;
dt=6e-12;         %time increment
T=10*2*L*n/c;     %total time

fmax = 1e9;       %maximun frequency
%fs=80*fmax;

TA=-T/2:dt:T/2;   %time axis for the signal
fs=1/dt;          %sampling frequency
%t = (-T/2/dt:1:T/2/dt)*dt;
Nt=round(T/dt);   

vsine = 1;
phi = vsine*sin(2*pi*fmax*TA);
EL1t=1.274e7*exp(1i*phi);                   %original signal
%plot(TA,(EL1t));
%FA = ((0:Nt-1)-floor(Nt/2))/Nt*fs;
FA = (-Nt/2:Nt/2-1)/Nt*fs;  %frequency axis
FP=fftshift(abs(fft(EL1t/Nt)));                %fft of the signal 
%fs=1/dt/Nt;
figure;
Z=plot(FA,FP);%fft plot with frequency axis

EL1tn=sqrt((T*P0/2/n/c/eps0)/(trapz(TA,abs(EL1t).^2)))*EL1t;  %original signal multiplied by normalization constant
Area=trapz(TA,2*n*c*eps0*abs(EL1tn).^2)/T                     %area under the curve(average power of the signal)

Area = 100.0000

I have an input of average power of 100W and I have normalized it as shown above.

But i want to do the same thing by converting the input signal(EL1t) into frequency axis(through fft) and then normalize it in frequency domain such that avearge power is 100W.

[Sam Chak]

Do you mean to integrate these 7 significant pillars (spectrums) over a certain frequency range?

P0=100;
eps0=8.854e-12;
L=10;
n=1.45;
Omega=1.020663297455751e+11;
GAMMAb=3.590391604102621e+08;
c=2.9979e8;
dt=6e-12;         %time increment
T=10*2*L*n/c;     %total time
fmax = 1e9;       %maximun frequency
TA=-T/2:dt:T/2;   %time axis for the signal
fs=1/dt;          %sampling frequency
Nt=round(T/dt);   
vsine = 1;
phi = vsine*sin(2*pi*fmax*TA);
EL1t=1.274e7*exp(1i*phi);                   %original signal
FA = (-Nt/2:Nt/2-1)/Nt*fs;  %frequency axis
FP=fftshift(abs(fft(EL1t/Nt)));                %fft of the signal 
figure;
Z=plot(FA,FP);%fft plot with frequency axis
xlim([-5e9 5e9]), grid on

When zooming in on the most prominent spectrum, a perfect triangular shape is seen. Shouldn't that resembles a Gaussian distribution curve or a Laplace distribution curve?

figure;
Z=plot(FA,FP);%fft plot with frequency axis
xlim([-2e6 2e6]), grid on

Answer 1

[Avni Agrawal]

Hi @Yogesh,

I understand that you are trying to normalize the average power of a signal to 100W in the frequency domain. To achieve this:

• Calculate the signal's Fourier Transform to represent it in the frequency domain.
• Determine the Power Spectral Density (PSD) by squaring the magnitude of the Fourier Transform and dividing by the number of points. Sum up the PSD across all frequencies to find the total power.
• If this total power isn't 100W, calculate a normalization factor that scales the total power to 100W. This factor is the square root of the desired power divided by the calculated total power. Multiply the time-domain signal by this normalization factor, and then recalculate its Fourier Transform and PSD.

The resulting spectrum should correspond to a signal with an average power of 100W in the frequency domain.

Please take a look at code snippet below:

clear all
close all
clc
P0=100;
eps0=8.854e-12;
L=10;
n=1.45;
Omega=1.020663297455751e+11;
GAMMAb=3.590391604102621e+08;
c=2.9979e8;
dt=6e-12;         %time increment
T=10*2*L*n/c;     %total time

fmax = 1e9;       %maximun frequency
%fs=80*fmax;

TA=-T/2:dt:T/2;   %time axis for the signal
fs=1/dt;          %sampling frequency
%t = (-T/2/dt:1:T/2/dt)*dt;
Nt=round(T/dt);   

vsine = 1;
phi = vsine*sin(2*pi*fmax*TA);
EL1t=1.274e7*exp(1i*phi);                   %original signal
%plot(TA,(EL1t));
%FA = ((0:Nt-1)-floor(Nt/2))/Nt*fs;
FA = (-Nt/2:Nt/2-1)/Nt*fs;  %frequency axis
FP=fftshift(abs(fft(EL1t/Nt)));                %fft of the signal 
%fs=1/dt/Nt;
figure;
Z=plot(FA,FP);%fft plot with frequency axis

EL1tn=sqrt((T*P0/2/n/c/eps0)/(trapz(TA,abs(EL1t).^2)))*EL1t;  %original signal multiplied by normalization constant
Area=trapz(TA,2*n*c*eps0*abs(EL1tn).^2)/T                     %area under the curve(average power of the signal)

Area = 100.0000

% Assuming EL1t, TA, and FA are defined as in your code snippet

% Step 2: Calculate PSD
PSD = (abs(fftshift(fft(EL1t))).^2)/Nt;

% Step 3: Calculate Total Power in Frequency Domain
% Since the signal is real, the power in the negative frequencies is mirrored in the positive frequencies.
% We can calculate the power in one side and multiply by 2 if needed, or just sum up the total power directly.
TotalPowerFreqDomain = sum(PSD)*(fs/Nt);  % fs/Nt is the frequency bin width

% Step 4: Normalize to Desired Average Power (100W)
DesiredPower = 100;  % Watts
NormalizationFactor = sqrt(DesiredPower/TotalPowerFreqDomain);
EL1t_normalized = EL1t * NormalizationFactor;

% Recalculate FFT and PSD for normalized signal
EL1t_normalized_FFT = fftshift(fft(EL1t_normalized));
PSD_normalized = (abs(EL1t_normalized_FFT).^2)/Nt;
NormalizedPowerFreqDomain = sum(PSD_normalized)*(fs/Nt);  % Should be close to 100W

% Plotting for verification
figure;
plot(FA, abs(EL1t_normalized_FFT)/Nt);  % Plot normalized spectrum
title('Normalized Spectrum');

% Check the normalized power
fprintf('Normalized Power in Frequency Domain: %f W\n', NormalizedPowerFreqDomain);

Normalized Power in Frequency Domain: 100.000000 W

I hope this helps!