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

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)

Answer 1

[Star Strider]

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.

.

[ime Rasenberg]

hey thanks for your answer, but I think I already have that, Ill show my code.

---

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)

---

if you run this code you will see one peak, that peak should consist of a few peaks, but i dont have the resolution in that erea to see what is going on with the harmonics, what could I do to get the FFT looking at the domain {f0-10a to f0+10a} so I can actualy see the peaks?

[Star Strider]

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.

.

Answer 2

[Walter Roberson]

https://www.mathworks.com/help/matlab/ref/double.nufft.html see the f parameter

This does not require uniform tspan since you can pass in the output t as the second parameter