Преобразование Фурье простой синусоидальной волны в Matlab - PullRequest
Новая
       26

Преобразование Фурье простой синусоидальной волны в Matlab

3 голосов
/ 02 апреля 2012

Я пытаюсь показать спектр простой синусоидальной волны, поскольку мы знаем, что одиночная синусоидальная волна с фиксированной частотой должна иметь пик в своем спектре. Я пишу этот код, но я не могу получить этот единственный пик, что неправильно в моем коде: 

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))

The result is here

Ответы [ 2 ]

2 голосов
/ 02 апреля 2012

Вы можете рассмотреть возможность использования встроенного преобразования Фурье, которое предоставляет MATLAB, вместо того, чтобы писать свое собственное. Смотри http://www.mathworks.se/help/techdoc/math/brentm1-1.html

Обновление

В функции fft есть несколько специфических особенностей, о которых вы должны знать. Во-первых, генерируемый массив необходимо нормализовать с коэффициентом 1/N, где N - это размер массива (это связано с реализацией функции MATLAB). Если вы этого не сделаете, все амплитуды в частотной области будут в N раз больше, чем они "на самом деле".

Во-вторых, вам нужно как-то найти частоты, которым соответствует каждый элемент в выходном массиве. Третья строка в приведенном ниже коде преобразует массив времен выборки в частоты, которым соответствует массив fourier.

Эта функция принимает сигнал во временной области и отображает его в частотной области. t - это временной массив, а y - это амплитуды сигнала. 

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');
0 голосов
/ 01 июня 2012

Ключевым моментом здесь является использование оконной функции. В результате код выглядит так: 

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))

И результат:

enter image description here

...