How can I define proper frequency domain of FFT function in Matlab

I'm having trouble of using function FFT and defining frequency domain

clear;
clc;
lambda_o=1e-6;    % 1[um] of wavelength
c=3e8;            % speed of light 
f_o=c/lambda_o;  % frequency (300 THz)
t=(-40:0.1:40)*1e-15;   % time domain ( -40fs ~ 40fs with dt=0.1fs)
dt=mean(diff(t));
tau=18*1e-15;           % Tau 18 fs
E_t=exp(-1i*2*pi*f_o*t).*exp(-pi*(t/tau).^2); % function that I'm interested
n=2^nextpow2(length(t)); %n=1024
Fs=800*1e12      %Sample Frequency (800THz)
df=Fs/n;         %Frequency resolution of 800THz/1024 (about 0.78125 THz)
E_f=fft(E_t,n);
E_f=E_f(n/2:n-1);  %Removing half-left elements of fourier transformed E_t
f=-Fs/2:df:Fs/2-df; %Definition of frequency domain (-400THz ~ 400THz)
f=f(n/2:n-1);  %Removing minus frequency domain
plot(f,abs(E_f)); %plotting the frequency spectrum

Thought the code above would make spike at 300THz frequency point since the function has a form of 300THz frequency.

But as you can see from the plot above it doesn't.

I suspect that this problem comes from the misunderstanding of defining frequency domain.

What am I missing? Please give me help.

And also Thank you for your kind comments in advance.

Answer 1

[Mohamad]

Hi , I correct the code , check

clear
clc
close all
lambda_o=1e-6;    % 1[um] of wavelength
c=3e8;            % speed of light 
f_o=c/lambda_o;  % frequency (300 THz)
Fs=25*f_o;      % Sampling Frequency 
Ts=1/Fs;        % Sampling Interval
femto=10^-15; % 1 Femto second
t=(-40*femto:Ts:40*femto);   % time domain ( -40fs ~ 40fs )
tau=18*femto;           % Tau 18 fs
E_t=exp(-1i*2*pi*f_o*t).*exp(-pi*(t/tau).^2); % function that I'm interested
n=2^nextpow2(length(t)); %
df=Fs/n;         % Frequency resolution of 800THz/1024 (about 0.78125 THz)
E_f=fft(E_t,n);
E_f=E_f(1:n/2+1);
E_f(2:end-1)=2*E_f(2:end-1);
f=(0:n/2)*Fs;
plot(f/1000,abs(E_f)); % plotting the frequency spectrum
grid ; xlabel(' frequency in GHz ') ; ylabel(' Magnitude ' )  

[재원 김]

Dear Mohamad, thanks for the comment

But unfortunately code didn't solve the problem.

I post the plot result of your code.

It doesn't show peak at the point of 300THz frequency.

Answer 2

[David Goodmanson]

Probably a bit late in the game, but the followiing code shows the peak at 300Thz. A couple of issues here: first to (a) get a frequency grid that shows the peak, and (b) if it is deemed to be worthwhile, to show the peak right at 300Thz rather than something fairly close.

I think the very first thing to do is NOT go to 2^n points with nextpow2.

All that does is confuse the situation, making it harder to accomplish (a), It also mucks up the resulting frequency grid, which is not commensurate with a frequency of 300THz, so you can't accomplish (b). The code below lops off the last point in the time domain to get an array of 600 points rather than 601. That allows a commensurate frequency grid, with one point at exactly 300THz. The fft is done on those 600 points.

The 2^n thing is long-established folklore and is done because of a supposed speed advantage. For a grid of 600 points, which is a highly composite number, there is no speed advantage at all to zeropadding 600 points to 1024 points before doing the fft. The supposed speed advantage is nonexistent. In fact, 1024 is slower. For 601 points, which is a prime number, yes, the fft on 1024 points is faster. But on my pc the fft for 601 points takes all of 20 microseconds, so any 2^n speed advantage is inconsequential. Avoiding ffts on arrays with prime length is important for large arrays, but not here.

clear
clc
close all
lambda_o=1e-6;    % 1[um] of wavelength
c=3e8;            % speed of light 
f0=c/lambda_o;  % frequency (300 THz)
Fs=25*f0;      % Sampling Frequency 
Ts=1/Fs;        % Sampling Interval
femto=10^-15; % 1 Femto second
t=(-40*femto:Ts:40*femto);   % time domain ( -40fs ~ 40fs )
% eliminate the last time point
t(end) = [];    
tau=18*femto;           % Tau 18 fs
E_t=exp(-1i*2*pi*f0*t).*exp(-pi*(t/tau).^2); % function that I'm interested
figure(1)
ploc(t,E_t); grid on
n = numel(t);
delf = Fs/n;         % freq grid spacing
% shift t = 0 to first point in time array;
% shift f = 0 to the middle of the frequency array
ampl = fftshift(fft(ifftshift(E_t)))/n; 
f = (-n/2:n/2-1)*delf;
figure(2)
plot(f,abs(ampl),'-o'); grid on
xlim([-800e12 200e12])