Trying to ubderstand the power distribution in fft plot
2 views (last 30 days)
Show older comments
clear all
close all
clc
L=10;
n=1.45;
c=2.9979e8;
dt = 6e-12;
T=10*2*L*n/c;
eps0=8.854e-12;
A=80e-12;
t = (-T/2/dt:1:T/2/dt)*dt;
Nt=round(T/dt);
fsine = 1e9;
vsine = 1;
phi = vsine*sin(2*pi*fsine*t);
EL1t=1.274e7*exp(1i*phi);
FP=fft(phi);
fs=1/dt/Nt;
Fs=(-1/dt/2:fs:1/dt/2-1);
figure
Z=plot(Fs,fftshift(abs(fft(EL1t/Nt).^2*2*n*c*eps0*A)));
As you see from the obtained fft plot , the peak of the graph at 0Hz is around 60W , but I am struggling to understand how the power is distributed throughout the plot.
The given input power is 100W I think. Shouldn't the central frequency at 0Hz be around 100W.
Where did the rest of power is what I am not understanding...
2 Comments
David Goodmanson
on 29 Apr 2024
Hello Yogesh,
the signal you are taking the fft of is proportional to
exp(i*const*sin(const*t))
is that your intent?
Accepted Answer
Mathieu NOE
on 30 Apr 2024
hello
the total power of the fft spectrum equals ~100W
pow_spectrum_total = 99.9502
but your fft spectrum spread it on multiple fft bins , and you are focusing only on the peak , not on the sum
L=10;
n=1.45;
c=2.9979e8;
dt = 6e-12;
T=10*2*L*n/c;
eps0=8.854e-12;
A=80e-12;
t = (-T/2/dt:1:T/2/dt)*dt;
Nt=round(T/dt);
fsine = 1e9;
vsine = 1;
phi = vsine*sin(2*pi*fsine*t);
EL1t=1.274e7*exp(1i*phi);
FP=fft(phi);
fs=1/dt/Nt;
Fs=(-1/dt/2:fs:1/dt/2-1);
pow_spectrum = fftshift(abs(fft(EL1t/Nt).^2*2*n*c*eps0*A));
pow_spectrum_total = sum(pow_spectrum)
figure
Z=plot(Fs,pow_spectrum);
More Answers (1)
David Goodmanson
on 30 Apr 2024
Edited: David Goodmanson
on 30 Apr 2024
Hello Yogesh,
The fact that all the abs(peaks)^2 add up to the expected total is just Parseval's theorem and has really nothing to do with the cause of the peaks. There are a couple of things going on here. The first is that the chosen frequency fsine = 1e9 does not have an integral number of cycles within the time window, and so it is not periodic in the time window.
t = 6e-12;
T=10*2*L*n/c;
t = (-T/2/dt:1:T/2/dt)*dt;
Nt =round(T/dt); % Nt = 161224
dt*Nt*fsine % cycles
ans = 967.3440
Just that situation all by itself means that a single peak in the frequency domain is not possible and there will be side peaks that mess up the result. By this I mean if you do an fft on
exp(2*pi*i*fsine*t) (1)
you will get a main peak and several smaller side peaks. This is because fsine = 1e9 is not commensurate with the frequencies in the fft. The code below replaces the time and frequency arrays with some similar ones for which 1e9 does fit exactly, having 1000 cycles. Then for the fft of (1) you get a single peak at the correct ampitude and no unwanted side peaks at all. After that change the fft calculation goes on just as before.
The more fundamental aspect is that you are transforming a function of the form
exp(vsine*sin(2*pi*fsine*t))
and its fourier transform definitely does have sidebands. The following expression gives their magnitude
exp(i*z*sin(theta)) = J_0(z)
+ Sum{k=1 to inf} J_2k(z)*[exp(i*2*k*theta) + exp(-i*2*k*theta)]
+ Sum{k=1 to inf} J_2k+1(z)*[exp(i*(2*k+1)*theta) - exp(-i*(2*k+1)*theta)]
Using
z = vsine and theta = 2*pi*fsine*t
you can identify the harmonics, and after collecting the all constants into C you can get analytic predictions for the peaks, which are shown with the red circles in the plot.
L=10;
n=1.45;
c=2.9979e8;
eps0=8.854e-12;
A=80e-12;
N = 2e5;
delt = 5e-12;
t = (-Nt/2:Nt/2-1)*delt;
delf = 1/(Nt*delt);
f = (-Nt/2:Nt/2-1)*delf;
fsine = 1e9; % N*delt*fsine = 1000 cycles
vsine = 1;
phi = vsine*sin(2*pi*fsine*t);
EL1t=1.274e7*exp(1i*phi);
figure(1)
plot(f,fftshift(abs(fft(EL1t/N).^2*2*n*c*eps0*A)));
hold on
C = 1.274e7^2*2*n*c*eps0*A
P0 = C*besselj(0,vsine)^2
P1 = C*besselj(1,vsine)^2
P2 = C*besselj(2,vsine)^2
fpks = [-2 -1 0 1 2]*1e9;
plot(fpks,[P2 P1 P0 P1 P2],'or')
xlim(1e11*[-.1 .1])
grid on
hold off
******************** also, note that
P0 +2*P1 +2*P2
ans = 99.8724
0 Comments
See Also
Categories
Find more on Fourier Analysis and Filtering in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!