how to define a frequenty domain to do an FFT over?

7 views (last 30 days)
At the moment the spectra I get from the FFT(x) function is in a too big domain.
I have one frequency at 5*10^14 and a second that is "a" (5*10^8) higher, they interact in a differential equation and generate harmonics, I am interested in those harmonics and how they evolve as I change "a". I was wondering how I can get a high-resolution FFT from f1 to f2 in high resolution so I can study this differential equation better. because now all the harmonics fall into one peak in the regular FFT.
Or do I have to program my own FFT to get it to be in this specific frequency domain.
my code looks as follows,
---clc, clear all, close all
%defining constants for differential equation
x0=-2;
x1=10^-12;
x2=10^-12;
%defining wavelengths
c=299792458;
w=550*10^-9;
f0=c/w;
a=500*10^6;
%defining time domain
te=10*x1;
dt=1*10^-16;
t=[0:dt:te];
%defining wave input
A1=0.06;
A2=0.2;
in=( A1 + A2 * exp( 1i*a*t ) ) .*exp( 1i*f0*t );
In=in.*in;
x=-1.9;
time=1;
%differential equation solved using eulers methode.
while te>t(time)
dx = ( ( x0-x(time) ) / x1 - In(time) ./ x2 * (exp( x(time) )-1) )*dt;
x(time+1)=x(time)+dx;
time=time+1;
end
out=in .* exp( ( 1 - 1i ) .* x ); %transforming wave equation
%FFT
L=length(out);
out2=out(L*3/4:1:L);
t2=t(L*3/4:1:L);
L=length(out2);
Y=fft(out2);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
P1=log(P1/max(P1));
f = 1/dt*(0:(L/2))/L;
plot(f,P1)

Accepted Answer

Star Strider
Star Strider on 16 Sep 2022
It would help to have an example.
The fft function requires that the original data be sampled at regular intervals. The easiest way to do that with the MATLAB ordinary differential equaiton integrators is to define ‘tspan’ as a vector of multiple time instants, rather than the beginning and ending times.
Example —
t_start = 0;
t_end = 25;
NrTimes = 50;
tspan = linspace(t_start, t_end, NrTimes);
If ‘NrTimes’ is a power of 2, the fft calculation will be faster, so:
NrTimes = 50;
NrTimes = 2^nextpow2(NrTimes);
tspan = linspace(t_start, t_end, NrTimes);
will also improve (increase) the fft frequency resolution.
.
  2 Comments
Star Strider
Star Strider on 16 Sep 2022
I am not certain what your code does, so it is difficult for me to make specific recommendations. (I thought you were using one of the MATLAB ODE integrators.)
One option is to increase the frequency resolution of the fft by zero-padding it. (Zero-padding it to a length that is an integer power-of-2 also increases its efficiency.) I did that here, however other than the expected ‘lobulation’ irregularities that result, I can see only one principal peak, even when using xlim to zoom in on the peak.
%defining constants for differential equation
x0=-2;
x1=10^-12;
x2=10^-12;
%defining wavelengths
c=299792458;
w=550*10^-9;
f0=c/w;
a=500*10^6;
%defining time domain
te=10*x1;
dt=1E-16;
t=[0:dt:te];
%defining wave input
A1=0.06;
A2=0.2;
in=( A1 + A2 * exp( 1i*a*t ) ) .*exp( 1i*f0*t );
In=in.*in;
x=-1.9;
time=1;
%differential equation solved using eulers methode.
while te>t(time)
dx = ( ( x0-x(time) ) / x1 - In(time) ./ x2 * (exp( x(time) )-1) )*dt;
x(time+1)=x(time)+dx;
time=time+1;
end
out=in .* exp( ( 1 - 1i ) .* x ); %transforming wave equation
%FFT
L=length(out);
out2=out(floor(L*3/4):1:L);
t2=t(floor(L*3/4):1:L);
figure
subplot(2,1,1)
plot(t2, real(out2))
grid
xlim([8 8.1]*1E-12) % Zoom
xlabel('Time')
ylabel('Amplitude')
title('Real')
subplot(2,1,2)
plot(t2, imag(out2))
grid
xlim([8 8.1]*1E-12) % Zoom
xlabel('Time')
ylabel('Amplitude')
title('Imaginary')
sgtitle('Time Domain Plot Of ‘out2’')
L=length(out2);
NFFT = 2^nextpow2(L);
Y=fft(out2,NFFT);
P2 = abs(Y/L);
P1 = P2(1:floor(L/2+1));
P1(2:end-1) = 2*P1(2:end-1);
P1=log(P1/max(P1));
% f = 1/dt*(0:(L/2))/L;
f = linspace(0, 1, numel(P1))/(2*dt);
figure
plot(f,P1)
xlim([1 1.3]*1E+14)
xlabel('Frequency')
ylabel('Magnitude')
The original signals appear to me to be regular sine and cosine curves with no harmonics or distortion.
I doubt any further increas in resolution is possible.
.

Sign in to comment.

More Answers (1)

Walter Roberson
Walter Roberson on 16 Sep 2022
Edited: Walter Roberson on 19 Sep 2022
This does not require uniform tspan since you can pass in the output t as the second parameter

Tags

Products


Release

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!