I try to show spectrum of simple sin wave as we know a single sin wave with fixed frequency must have on peak in its spectrum I write this code but I can't get this one peak what is wrong in my code:
clc
nsteps=200;%number of signal elements in time domain
i=sqrt(-1);
NFREQS=100;%number of elements in frequency domain
ddx=1e-9;
dt=ddx/(6e8);%separation between each time domain elements
lambdai=150e-9;
lambdaf=500e-9;
freqi=3e8/lambdai;
freqf=3e8/lambdaf;
freq=zeros(1,NFREQS);
for j=1:NFREQS
freq(j)=freqi-j*(freqi-freqf)/NFREQS;%desired frequency domain
end
arg=2*pi*freq*dt;
et=zeros(nsteps,1);
for j=1:nsteps
et(j)=sin(2*pi*3e15*j*dt);%sin wave in time domain
end
e=zeros(NFREQS,1);
for n=1:NFREQS
for j=1:nsteps
e(n)=e(n)+et(j)*exp(-i*arg(n)*n);%sin wave in frequency domain
end
end
lambda=linspace(lambdai,lambdaf,NFREQS);
plot(lambda,abs(e))
You could consider using the built in Fourier transform that MATLAB provides instead of writing your own. See http://www.mathworks.se/help/techdoc/math/brentm1-1.html
Update
There are a few peculiar things about the fft
function you should know about. The first thing is that the array that it generates needs to be normalized by a factor of 1/N
where N
is the size of the array (this is because of the implementation of the MATLAB function). If you don't do that all the amplitudes in frequency domain will be N
times larger than they "actually" are.
The second thing is that you need to somehow find the frequencies that each element in the output array corresponds to. The third line in the code below converts the array of sample times to the frequencies that the fourier
array corresponds to.
This function takes a signal in time domain and plots it in frequency domain. t
is the time array and y
is the signal amplitudes.
function plotSpectrum(t, y)
freqUnit = 1 / (abs(t(2) - t(1)) * length(t));
f = -length(t) * freqUnit : freqUnit : (length(t) - 1) * freqUnit;
oneSidedFFT = fft(y);
fourier = horzcat(oneSidedFFT, oneSidedFFT);
plot(f, abs(fourier));
xlabel('Frequency');
ylabel('Magnitude');
The key point in here is using window function. The resulted code is this:
clc
clear
nsteps=40000;
i=sqrt(-1);
Sc=0.5;
ddx=1e-9;
dt=ddx*Sc/(3e8);
lambdai=100e-9;
lambdaf=700e-9;
lambda=lambdai:1e-9:lambdaf;
[m,NFREQS]=size(lambda);
freq=3e8./lambda;
et=zeros(nsteps,1);
for j=1:nsteps
et(j)=sin(2*pi*3e8/(300e-9)*dt*j);%sin wave in time domain
end
w=blackman(nsteps);%window function for periodic functions
for j=1:nsteps
et(j)=et(j)*w(j);
end
e=zeros(1,NFREQS);
for n=1:NFREQS
for j=1:nsteps
e(n)=e(n)+et(j)*exp(-i*2*pi*freq(n)*dt*j);
end
end
plot(lambda,abs(e))
And the result is: