頻率抽取法FFT的程序

clear all;

%輸入f、N、T、是否補零（補幾個零）

f=input('Input frequency of the signal: f\n');

N=input('Input number of pointsl: N\n');

T=input('Input sampling time: T\n');

flag=input('Add zero too sampling signal or not? yes=1 no=0\n');

if(flag)

ZeroNum=input('Input nmber of zeros\n');

else

ZeroNum=0;

end

%生成信號，signal是原信號。signal為采樣信號。

fs=1/T;

t=0:0.00001:T*(N+ZeroNum-1);

signal=sin(2*pi*f*t);

t2=0:T:T*(N+ZeroNum-1);

signal2=sin(2*pi*f*t2);

if (flag)

signal2=[signal2 zeros(1, ZeroNum)];

end

%畫出原信號及采樣信號。

figure;

subplot(2,1,1);

plot(t,signal);

xlabel('Time(s)');

ylabel('Amplitude(volt)');

title('Singnal');

hold on;

subplot(2,1,1);

stem(t2,signal2,'r');

axis([0 T*(N+ZeroNum) -1 1]);

%作FFT變換，計算其幅值，歸壹化處理，並畫出頻譜。

Y = fft(signal2,N);

Pyy = Y.* conj(Y) ;

Pyy=(Pyy/sum(Pyy))*2;

f=0:fs/(N-1):fs/2;4

subplot(2,1,2);

bar(f,Pyy(1:N/2));

xlabel('Frequency(Hz)');

ylabel('Amplitude');

title('Frequency compnents of signal');

axis([0 fs/2 0 ceil(max(Pyy))])

grid on;