Frenquency shift between the value and plot

Hello,

I wrote a programm that is supposed to do the FFT of a signal. The problem I have is that the peak is of the curve should line up with the orange one (i will send the program). But in fact they do not line up and I do not know why. There is a small shift between the value of the frequency 2.849*10^5GHz and the value displayed in plot (2.848*10^5 GHz). I do not know the origin of this small shift.

Thanks for your help

Answer 1

[Paul]

Hi Achraf,

The signal p is intended to have frequency content at f1 and f2.

lbd1  = 1053;
f1    = 3e8/(lbd1*10^-9)

f1 = 2.8490e+14

and f2 = f1.

However the sampling frequency is

sample_rate = 1e12;

which is two orders of magnitude smaller than f1 and f2, when it should be more than 2x larger.

Consequently, the frequency in p aliases down to

aliasfreq =  9.9715e+10;

which is the offset in the plot (before the divide by 10e8)

f1 - aliasfreq

ans = 2.8480e+14

TBH, I had trouble following the code and all of the offsets and manipulations, but I'm pretty sure this aliasing is a problem.

Also, because N is odd for this problem, the f vector should be

f = (-((N-1)/2):((N-1)/2))/N*sample_rate + f0;

which will be slightly different than f in the code.

[Achraf AYEB]

Hello,

I have tried the alliasing problem but it doesnot work. I still have the same problem. I don't know if it is some fft that is not working as it should be. When I increase the sampling frequency it is even worse here it is 1e13

here it is 1e12

[Paul]

Can you post the updated code? Please not as an .mlx file. Either attach as an m-file, or copy/paste into a comment in this thread.

[Achraf AYEB]

clear all

close all

sample_rate=1e13;%Sample rate

T=1/sample_rate;

L=1e6; %signal length

t =(0:L)*T;%time vector

f0=3e8/1053e-9;

%laser frequency Hz => lambda=C/F=C.T

w=2*pi*f0;%laser pulsation

fm=3e9;%modulation frequency in GHz

wmod=2*pi*fm;%phase modulation

N=length(t);

lbd1=1053;

f1=3e8/(lbd1*10^-9); w1=2*pi*f1;

lbd2=1053;

f2=3e8/(lbd2*10^-9); w2=2*pi*f2;

m=0;%modulation depth

gaussian_order=0;

g0=6.6;%ampli gain

l=2; %material length

shift=40e-9;%temporal shift

shift2=80e-9;%temporal shift

sigma=10e-9;%Gaussian HWHM

p=exp(-1/2*((t-shift)/(sigma/2)).^gaussian_order).*exp(1i.*w1.*t)+exp(-1/2*((t-shift2)/(sigma/2)).^gaussian_order).*exp(1i.*w2.*t); %square pulse modelised by a super-gaussian

x=p.*exp(1i*m*cos(wmod*t));

y=fft(x);%fft supergaussian

y=y/max(abs(y));%normalisation fft

z=fftshift(y);% rearranges a Fourier transform X by shifting the zero-frequency component to the center of the array

f = (-((N-1)/2):((N-1)/2))/N*sample_rate + f0;

%Frequency vector

n=length(f);

%number of element in f vector

%Material fluorescence

fluo=g0*exp(-1/2*((3e8./f-1053e-9)/0.7e-9).^2);

amp=sqrt(abs(z));

%norm of z

gain=amp.*exp(l*fluo);

%total gain

gain=gain/max(gain);

yamp=fftshift(z.*exp(l*fluo));% rearranges a Fourier transform X by shifting the zero-frequency component to the center of the array

p_freq=fft(p);%fft of square pulse

p_freqshifft=fftshift(p_freq);% rearranges a Fourier transform X by shifting the zero-frequency component to the center of the array

p_freqshifft=p_freqshifft/max(abs(p_freqshifft));%normalisation

%Time in ns

t=t*10e8;

%Frequency in GHz

f=f/10e8;

%Display graphs in 2*2

subplot(2,1,1)

%Display fluorescence

plot(f,fluo/max(fluo))

%normalized fluorescence

legend('Fluo')

hold

%impulsion spectrum

plot(f(1:n),abs(p_freqshifft(1:n)))

%normalized fluorescence

legend('Pulse without modulation')

hold

title('Power spectrum')

xlabel('Frequence (GHz)');

ylabel('normalized amplitude (U.A.)')

legend('Fluo')

%hold all

subplot(2,1,2)

plot(f(1:n),gain(1:n))

title('Power spectrum')

xlabel('Frequency (GHz)')

ylabel('Normalized Intensity (A.U.)')

[Paul]

In the current code, f1 = f2 = 2.849e14. So sample_rate needs to be larger than 2.849e14 by at least a factor of two. sample_rate of 1e12 and 1e13 are still too small. Once this issue is addressed, there's something else to consider.

Here's the code, but I made sample_rate approximately 5x larger than f1 and f2

%sample_rate=1e13;%Sample rate 
sample_rate=1.4e15;
T=1/sample_rate;
L=1e6; %signal length
t =(0:L)*T;%time vector
f0=3e8/1053e-9; 
%laser frequency Hz => lambda=C/F=C.T
%w=2*pi*f0;%laser pulsation
fm=3e9;%modulation frequency in GHz
wmod=2*pi*fm;%phase modulation
N=length(t); 
lbd1=1053;
f1=3e8/(lbd1*10^-9); w1=2*pi*f1;
lbd2=1053; 
f2=3e8/(lbd2*10^-9); w2=2*pi*f2;
m=0;%modulation depth
gaussian_order=0;
g0=6.6;%ampli gain
l=2; %material length
shift=40e-9;%temporal shift
shift2=80e-9;%temporal shift
sigma=10e-9;%Gaussian HWHM
p=exp(-1/2*((t-shift)/(sigma/2)).^gaussian_order).*exp(1i.*w1.*t)+exp(-1/2*((t-shift2)/(sigma/2)).^gaussian_order).*exp(1i.*w2.*t); %square pulse modelised by a super-gaussian
x=p.*exp(1i*m*cos(wmod*t)); 
y=fft(x);%fft supergaussian
y=y/max(abs(y));%normalisation fft
z=fftshift(y);% rearranges a Fourier transform X by shifting the zero-frequency component to the center of the array
f = (-((N-1)/2):((N-1)/2))/N*sample_rate + f0;
%Frequency vector
n=length(f); 
%number of element in f vector
%Material fluorescence
fluo=g0*exp(-1/2*((3e8./f-1053e-9)/0.7e-9).^2);
amp=sqrt(abs(z)); 
%norm of z
gain=amp.*exp(l*fluo); 
%total gain
gain=gain/max(gain);
%yamp=fftshift(z.*exp(l*fluo));% rearranges a Fourier transform X by shifting the zero-frequency component to the center of the array
p_freq=fft(p);%fft of square pulse
p_freqshifft=fftshift(p_freq);% rearranges a Fourier transform X by shifting the zero-frequency component to the center of the array
p_freqshifft=p_freqshifft/max(abs(p_freqshifft));%normalisation

Based on the defnition of p, we should expect its DFT to have a peak at 2.849e14 Hz, based on the values of f1 and f2.

freqp = (-((N-1)/2):((N-1)/2))/N*sample_rate;
figure
plot(freqp,abs(p_freqshifft)),xlabel('Hz')
xlim([2.8 2.9]*1e14)

The peak is exactly where it should be when plotting against fp, which is the correct frequency vector associated with p_freqshifft.

Let's keep everything in Hz

%Time in ns
%t=t*10e8; 

%Frequency in GHz
%f=f/10e8;

Now make the plots

figure
%Display graphs in 2*2
subplot(2,1,1)
%Display fluorescence
plot(f,fluo/max(fluo)) 
%normalized fluorescence
legend('Fluo')
hold

Current plot held

%impulsion spectrum
plot(f(1:n),abs(p_freqshifft(1:n))) 
%normalized fluorescence
legend('Pulse without modulation')
hold

Current plot released

title('Power spectrum')
xlabel('Frequence (GHz)');
ylabel('normalized amplitude (U.A.)')
legend('Fluo','p_freqshifft','Interpreter','none')
%hold all
subplot(2,1,2)
plot(f(1:n),gain(1:n))
title('Power spectrum')
xlabel('Frequency (GHz)')
ylabel('Normalized Intensity (A.U.)')

Now, the peak of p_freqshifft appears to be at the wrong location. But that happens because because it's being plotted against f, which is shifted by f0 (which, coincidentally or not is equal to f1 and f2) from freqp, so the peak appears at

fpeak = f0 + f1

fpeak = 5.6980e+14

figure
plot(f(1:n),abs(p_freqshifft(1:n)))
xline(fpeak,'r')
xlim([5.695 5.705]*1e14)

[Achraf AYEB]

Thank you Paul,

I understand now, in a way I have to make an x-axis for each one the curve. Have a nice day.