f = 1/(200e-9);? ?%方波信號(hào)的頻率，5Mhz，200ns

L = 512;? ? ? ? ? ? ?% 每個(gè)周期采樣的點(diǎn)數(shù)

Fs = L*f;? ? ? ? ? ? % Sampling frequency ，采樣率=信號(hào)頻率*每個(gè)周期的采樣數(shù)? ? ? ? ? ? ? ? ? ?

T = 1/Fs;? ? ? ? ? ? ?% Sampling period? ? ?兩次采樣的時(shí)間間隔

t = (0:L-1)*T;? ? ? ? % Time vector? ? 通過采樣時(shí)間間隔和采樣的點(diǎn)數(shù)生成時(shí)間矩陣。這里選取一個(gè)周期的方波信號(hào)：512個(gè)點(diǎn)=1個(gè)方波周期，時(shí)間從0開始所以最后減去一個(gè)點(diǎn)。

x = square(t*f*2*pi); %square可以方便的生成方波信號(hào)，生成方式和正弦信號(hào)一樣，角速度*時(shí)間，官方說明文檔有很多實(shí)例。

plot(t,x);? %畫出方波信號(hào)幅值隨時(shí)間的變化圖

%n = 2^nextpow2(L); %在有些文檔里對(duì)進(jìn)行快速傅里葉變化的點(diǎn)數(shù)進(jìn)行補(bǔ)零。（可能可以加快運(yùn)算，測試了1024和1000點(diǎn)精度上沒有區(qū)別）

Y = fft(x); %快速傅里葉變換

P2 = abs(Y/L); %雙邊頻譜，Matlab做快速傅里葉變換后會(huì)保留負(fù)頻部分需要后處理。這里取絕對(duì)值和除上數(shù)據(jù)長度，是為了還原縱坐標(biāo)和功率的對(duì)應(yīng)關(guān)系。（具體原有可能和傅里葉快速變換的數(shù)學(xué)推導(dǎo)有關(guān)）

P1 = P2(1:L/2+1); %截取頻譜一半加一個(gè)點(diǎn)（方法參考的官方文檔加一個(gè)點(diǎn)可能是由于直流部分和最高頻只有一個(gè)？）

P1(2:end-1) = 2*P1(2:end-1); %單邊頻譜，除頭尾乘2。根據(jù)上面猜想對(duì)有正負(fù)頻兩部分的進(jìn)行合并，還原正確的頻間功率關(guān)系

f_fft = Fs*(0:(L/2))/L;?%給這截取一半的頻譜圖對(duì)應(yīng)正確的橫坐標(biāo)，頻率。257個(gè)點(diǎn)有256個(gè)間隔。

plot(f_fft,P1) %畫出頻譜圖