fft演算法的matlab實現
⑴ matlab中fft（）函數是什麼意思
FFT（快速傅里葉變換）是一種實現DFT（離散傅里葉變換）的快速演算法，是利用復數形式的離散傅里葉變換來計算實數形式的離散傅里葉變換，matlab中的fft()函數是實現該演算法的實現。
MATLAB它將數值分析、矩陣計算、科學數據可視化以及非線性動態系統的建模和模擬等諸多強大功能集成在一個易於使用的視窗環境中，為科學研究、工程設計以及必須進行有效數值計算的眾多科學領域提供了一種全面的解決方案，並在很大程度上擺脫了傳統非互動式程序設計語言（如C、Fortran）的編輯模式，代表了當今國際科學計算軟體的先進水平。
快速傅里葉變換, 即利用計算機計算離散傅里葉變換（DFT)的高效、快速計算方法的統稱，簡稱FFT。快速傅里葉變換是1965年由J.W.庫利和T.W.圖基提出的。採用這種演算法能使計算機計算離散傅里葉變換所需要的乘法次數大為減少，特別是被變換的抽樣點數N越多，FFT演算法計算量的節省就越顯著。
(1)fft演算法的matlab實現擴展閱讀：
matlab優勢特點：
1、高效的數值計算及符號計算功能，能使用戶從繁雜的數學運算分析中解脫出來；
2、具有完備的圖形處理功能，實現計算結果和編程的可視化；
3、友好的用戶界面及接近數學表達式的自然化語言，使學者易於學習和掌握；
4、功能豐富的應用工具箱(如信號處理工具箱、通信工具箱等) ，為用戶提供了大量方便實用的處理工具。
參考資料來源：
網路-快速傅里葉變換
網路-MATLAB
⑵ fft演算法matlab的實現代碼！完整版的！
function result = MyFFT(vector)
result = fft(vector);
⑶ 如何應用matlab進行fft分析
FFT是離散傅立葉變換的快速演算法，可以將一個信號變換
到頻域。有些信號在時域上是很難看出什麼特徵的，但是如
果變換到頻域之後，就很容易看出特徵了。這就是很多信號
分析採用FFT變換的原因。另外，FFT可以將一個信號的頻譜
提取出來，這在頻譜分析方面也是經常用的。
雖然很多人都知道FFT是什麼，可以用來做什麼，怎麼去
做，但是卻不知道FFT之後的結果是什意思、如何決定要使用
多少點來做FFT。
現在圈圈就根據實際經驗來說說FFT結果的具體物理意義。
一個模擬信號，經過ADC采樣之後，就變成了數字信號。采樣
定理告訴我們，采樣頻率要大於信號頻率的兩倍，這些我就
不在此羅嗦了。
采樣得到的數字信號，就可以做FFT變換了。N個采樣點，
經過FFT之後，就可以得到N個點的FFT結果。為了方便進行FFT
運算，通常N取2的整數次方。
假設采樣頻率為Fs，信號頻率F，采樣點數為N。那麼FFT
之後結果就是一個為N點的復數。每一個點就對應著一個頻率
點。這個點的模值，就是該頻率值下的幅度特性。具體跟原始
信號的幅度有什麼關系呢？假設原始信號的峰值為A，那麼FFT
的結果的每個點（除了第一個點直流分量之外）的模值就是A
的N/2倍。而第一個點就是直流分量，它的模值就是直流分量
的N倍。而每個點的相位呢，就是在該頻率下的信號的相位。
第一個點表示直流分量（即0Hz），而最後一個點N的再下一個
點（實際上這個點是不存在的，這里是假設的第N+1個點，也
可以看做是將第一個點分做兩半分，另一半移到最後）則表示
采樣頻率Fs，這中間被N-1個點平均分成N等份，每個點的頻率
依次增加。例如某點n所表示的頻率為：Fn=(n-1)*Fs/N。
由上面的公式可以看出，Fn所能分辨到頻率為為Fs/N，如果
采樣頻率Fs為1024Hz，采樣點數為1024點，則可以分辨到1Hz。
1024Hz的采樣率采樣1024點，剛好是1秒，也就是說，采樣1秒
時間的信號並做FFT，則結果可以分析到1Hz，如果采樣2秒時
間的信號並做FFT，則結果可以分析到0.5Hz。如果要提高頻率
分辨力，則必須增加采樣點數，也即采樣時間。頻率解析度和
采樣時間是倒數關系。
假設FFT之後某點n用復數a+bi表示，那麼這個復數的模就是
An=根號a*a+b*b，相位就是Pn=atan2(b,a)。根據以上的結果，
就可以計算出n點（n≠1，且n<=N/2）對應的信號的表達式為：
An/(N/2)*cos(2*pi*Fn*t+Pn)，即2*An/N*cos(2*pi*Fn*t+Pn)。
對於n=1點的信號，是直流分量，幅度即為A1/N。
由於FFT結果的對稱性，通常我們只使用前半部分的結果，
即小於采樣頻率一半的結果。
好了，說了半天，看著公式也暈，下面圈圈以一個實際的
信號來做說明。
假設我們有一個信號，它含有2V的直流分量，頻率為50Hz、
相位為-30度、幅度為3V的交流信號，以及一個頻率為75Hz、
相位為90度、幅度為1.5V的交流信號。用數學表達式就是如下：
S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180)
式中cos參數為弧度，所以-30度和90度要分別換算成弧度。
我們以256Hz的采樣率對這個信號進行采樣，總共采樣256點。
按照我們上面的分析，Fn=(n-1)*Fs/N，我們可以知道，每兩個
點之間的間距就是1Hz，第n個點的頻率就是n-1。我們的信號
有3個頻率：0Hz、50Hz、75Hz，應該分別在第1個點、第51個點、
第76個點上出現峰值，其它各點應該接近0。實際情況如何呢？
我們來看看FFT的結果的模值如圖所示。
圖1 FFT結果
從圖中我們可以看到，在第1點、第51點、和第76點附近有
比較大的值。我們分別將這三個點附近的數據拿上來細看：
1點： 512+0i
2點： -2.6195E-14 - 1.4162E-13i
3點： -2.8586E-14 - 1.1898E-13i
50點：-6.2076E-13 - 2.1713E-12i
51點：332.55 - 192i
52點：-1.6707E-12 - 1.5241E-12i
75點：-2.2199E-13 -1.0076E-12i
76點：3.4315E-12 + 192i
77點：-3.0263E-14 +7.5609E-13i
很明顯，1點、51點、76點的值都比較大，它附近的點值
都很小，可以認為是0，即在那些頻率點上的信號幅度為0。
接著，我們來計算各點的幅度值。分別計算這三個點的模值，
結果如下：
1點： 512
51點：384
76點：192
按照公式，可以計算出直流分量為：512/N=512/256=2；
50Hz信號的幅度為：384/(N/2)=384/(256/2)=3；75Hz信號的
幅度為192/(N/2)=192/(256/2)=1.5。可見，從頻譜分析出來
的幅度是正確的。
然後再來計算相位信息。直流信號沒有相位可言，不用管
它。先計算50Hz信號的相位，atan2(-192, 332.55)=-0.5236,
結果是弧度，換算為角度就是180*(-0.5236)/pi=-30.0001。再
計算75Hz信號的相位，atan2(192, 3.4315E-12)=1.5708弧度，
換算成角度就是180*1.5708/pi=90.0002。可見，相位也是對的。
根據FFT結果以及上面的分析計算，我們就可以寫出信號的表達
式了，它就是我們開始提供的信號。
總結：假設采樣頻率為Fs，采樣點數為N，做FFT之後，某
一點n（n從1開始）表示的頻率為：Fn=(n-1)*Fs/N；該點的模值
除以N/2就是對應該頻率下的信號的幅度（對於直流信號是除以
N）；該點的相位即是對應該頻率下的信號的相位。相位的計算
可用函數atan2(b,a)計算。atan2(b,a)是求坐標為(a,b)點的角
度值，范圍從-pi到pi。要精確到xHz，則需要采樣長度為1/x秒
的信號，並做FFT。要提高頻率解析度，就需要增加采樣點數，
這在一些實際的應用中是不現實的，需要在較短的時間內完成
分析。解決這個問題的方法有頻率細分法，比較簡單的方法是
采樣比較短時間的信號，然後在後面補充一定數量的0，使其長度
達到需要的點數，再做FFT，這在一定程度上能夠提高頻率分辨力。
具體的頻率細分法可參考相關文獻。
[附錄：本測試數據使用的matlab程序]
close all; %先關閉所有圖片
Adc=2; %直流分量幅度
A1=3; %頻率F1信號的幅度
A2=1.5; %頻率F2信號的幅度
F1=50; %信號1頻率(Hz)
F2=75; %信號2頻率(Hz)
Fs=256; %采樣頻率(Hz)
P1=-30; %信號1相位(度)
P2=90; %信號相位(度)
N=256; %采樣點數
t=[0:1/Fs:N/Fs]; %采樣時刻
%信號
S=Adc+A1*cos(2*pi*F1*t+pi*P1/180)+A2*cos(2*pi*F2*t+pi*P2/180);
%顯示原始信號
plot(S);
title('原始信號');
figure;
Y = fft(S,N); %做FFT變換
Ayy = (abs(Y)); %取模
plot(Ayy(1:N)); %顯示原始的FFT模值結果
title('FFT 模值');
figure;
Ayy=Ayy/(N/2); %換算成實際的幅度
Ayy(1)=Ayy(1)/2;
F=([1:N]-1)*Fs/N; %換算成實際的頻率值
plot(F(1:N/2),Ayy(1:N/2)); %顯示換算後的FFT模值結果
title('幅度-頻率曲線圖');
figure;
Pyy=[1:N/2];
for i="1:N/2"
Pyy(i)=phase(Y(i)); %計算相位
Pyy(i)=Pyy(i)*180/pi; %換算為角度
end;
plot(F(1:N/2),Pyy(1:N/2)); %顯示相點陣圖
title('相位-頻率曲線圖');
看完這個你就明白諧波分析了。
⑷ FFT測量相位具體演算法。在matlab中如何使用進行編程
% 下面的程序里Pn 存的就是基波相位 如果求的是諧波相位，稍微修改即可
x = load('data.dat'); %load 數據
fs=10000; % 采樣頻率，自己根據實際情況設置
N=length(x); % x 是待分析的數據
n=1:N;
%1-FFT
X=fft(x); % FFT
X=X(1:N/2);
Xabs=abs(X);
Xabs(1) = 0; %直流分量置0
[Amax,index]=max(Xabs);
if(Xabs(index-1) > Xabs(index+1))
a1 = Xabs(index-1) / Xabs(index);
r1 = 1/(1+a1);
k01 = index -1;
else
a1 = Xabs(index) / Xabs(index+1);
r1 = 1/(1+a1);
k01 = index;
end
Fn = (k01+r1-1)*fs/N; %基波頻率
An = 2*pi*r1*Xabs(k01)/(N*sin(r1*pi)); %基波幅值
Pn = phase(X(k01))-pi*r1; %基波相位 單位弧度
Pn = mod(Pn(1),pi);
⑸ 用matlab如何實現fft變換
Matlab中FFT有1D和2D的，FFT得到的是信號的頻譜即t－》f
如
clear
%編寫駱遙
fs=1000
t=0:1/fs:0.6;
f1=100;
f2=300;
x=sin(2*pi*f1*t)+sin(2*pi*f2*t);
subplot(711)
plot(x);
title('f1（100Hz）\f2(300Hz)的正弦信號，初相0')
xlabel('序列（n）')
grid on
number=512
y=fft(x,number);
n=0:length(y)-1;
f=fs*n/length(y);
subplot(713)
plot(f,abs(y));
title('f1\f2的正弦信號的FFT（512點）')
xlabel('頻率Hz')
grid on
x=x+randn(1,length(x));
subplot(715)
plot(x);
title('原f1\f2的正弦信號（含隨機雜訊）')
xlabel('序列（n）')
grid on
y=fft(x,number);
n=0:length(y)-1;
f=fs*n/length(y);
subplot(717)
plot(f,abs(y));
title('原f1\f2的正弦信號（含隨機雜訊）的FFT（512點）')
xlabel('頻率Hz')
grid on
⑹ MATLAB中的FFT的采樣頻率和采樣點怎樣確定
在MATLAB中做FFT，首先編寫函數，對不同的采樣頻率和采樣點數，計算FFT後的頻率序列及其對應的幅值：
function[f amplitude]=yopheeFFT(sampleRate,FFT_points)
n=0:FFT_points-1;
t=n/sampleRate;%采樣時間序列
f_All=n*sampleRate/FFT_points;%頻率序列 %構造混有雜訊的周期信號並采樣
signal=2*sin(2*pi*10*t)+1*sin(2*pi*20.25*t)+0.2*randn(size(t));%對信號進行快速Fourier變換，並求振幅
amplitude_All=abs(fft(signal,FFT_points))*2/FFT_points;
f=f_All(1:FFT_points/2);
amplitude=amplitude_All(1:FFT_points/2);
(6)fft演算法的matlab實現擴展閱讀
MATLAB中FFT函數的意義：
FFT是離散傅立葉變換的快速演算法，可以將一個信號變換到頻域。有些信號在時域上是很難看出什麼特徵的，但是如果變換到頻域之後，就很容易看出特徵了。這就是很多信號分析採用FFT變換的原因。另外，FFT可以將一個信號的頻譜提取出來，這在頻譜分析方面也是經常用的。
模擬信號經過ADC采樣之後變成數字信號,可對此數字信號做FFT變換。N個采樣點經過FFT之後就可以得到N個點的FFT結果。為了方便進行FFT運算，通常N取2的整數次冪。
假設采樣頻率為Fs，信號頻率為F，采樣點數為N。則FFT之後結果為N點復數，其中每一個點對應著一個頻率點，該點復數的模值為原始信號在該頻率值下的幅度特性。
具體為：假設原始信號在某頻率點的幅值為A，則該頻點對應的FFT點復數的模值為A的N/2倍。而FFT第一點為原始信號的直流分量，其模值為原始信號模值的N倍。對於相位，FFT復數的相位即為原始信號在該頻率點處的相位。
⑺ 基-2fft演算法的軟體實現 matlab代碼
% 基於Matlab的時間抽取基2FFT演算法
function y=myditfft(x)
%本程序對輸入序列實現DIT-FFT基2演算法，點數取大於等於長度的2的冪次
%------------------------------------
% Leo's fft program(改編網上的一個程序)
%------------------------------------
m=log2(2^nextpow2(length(x))); %求的x長度對應的2的最低冪次m
N=2^m;
if length(x)<N
x=[x,zeros(1,N-length(x))]; %若長度不是2的冪，補0到2的整數冪
end
x;
%--------------------------------------------------------------------------
%對輸入序列進行倒序
%如果輸入序列的自然順序號I用二進制數（例如n2n1n0）表示
%則其倒位序J對應的二進制數就是（n0n1n2），這樣，在原來自然順序時應該放x(I)的
%單元，現在倒位序後應放x(J)。
%--------------------------------------------------------------------------
%以下程序相當於以下程序:
%nxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1; %求1：2^m數列的倒序
%y=x(nxd); %將倒序排列作為初始值
%--------------------------------------------------------------------------
NV2=N/2;
NM1=N-1;
I=0;
J=0;
while I<NM1
if I<J
T=x(J+1);
x(J+1)=x(I+1);
x(I+1)=T;
end
K=NV2;
while K<=J
J=J-K;
K=K/2;
end
J=J+K;
I=I+1;
end
x;
%--------------------------------------------------------------------------
%以下程序解釋:
%第一級從x(0)開始,跨接一階蝶形,再取每條對稱
%第二級從x(0)開始,跨接兩階蝶形,再取每條對稱
%第m級從x(0)開始,跨接2^(m-1)階蝶形,再取每條對稱....
%--------------------------------------------------------------------------
for mm=1:m %將DFT做m次基2分解，從左到右，對每次分解作DFT運算
Nmr=2^mm;
u=1; %旋轉因子u初始化
WN=exp(-j*2*pi/Nmr); %本次分解的基本DFT因子WN＝exp(-i*2*pi/Nmr)
for n=1:Nmr/2 %本次跨越間隔內的各次碟形運算
for k=n:Nmr:N %本次碟形運算的跨越間隔為Nmr=2^mm
kp=k+Nmr/2; %確定碟形運算的對應單元下標(對稱性)
t=x(kp)*u; %碟形運算的乘積項
x(kp)=x(k)-t; %碟形運算的加法項
x(k)=x(k)+t;
end
u=u*WN; %修改旋轉因子，多乘一個基本DFT因子WN
end
end
y=x; %輸出
⑻ 用matlab編寫實現fft的程序。
function y=myditfft(x)
%本程序對輸入序列實現DIT-FFT基2演算法，點數取大於等於長度的2的冪次
%------------------------------------
% myditfft.c
%------------------------------------
m=nextpow2(x); %求的x長度對應的2的最低冪次m
N=2^m;
if length(x)<N
x=[x,zeros(1,N-length(x))]; %若的長度不是2的冪，補0到2的整數冪
end
nxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1; %求1：2^m數列的倒序
y=x(nxd); %將倒序排列作為的初始值
for mm=1:m %將DFT做m次基2分解，從左到右，對每次分解作DFT運算
Nmr=2^mm;
u=1; %旋轉因子u初始化
WN=exp(-i*2*pi/Nmr); %本次分解的基本DFT因子WN＝exp(-i*2*pi/Nmr)
for j=1:Nmr/2 %本次跨越間隔內的各次碟形運算
for k=j:Nmr:N %本次碟形運算的跨越間隔為Nmr=2^mm
kp=k+Nmr/2; %確定碟形運算的對應單元下標
t=y(kp)*u; %碟形運算的乘積項
y(kp)=y(k)-t; %碟形運算的加法項
y(k)=y(k)+t;
end
u=u*WN; %修改旋轉因子，多乘一個基本DFT因子WN
end
end
⑼ 基2—fft演算法的軟體實現(MATLAB代碼)
參考網路： clc; clear all; close all; x=ones(1,128); %輸入的信號，自己可以改變 %整體運用原位計算 m=nextpow2(x);N=2^m; % 求x的長度對應的2的最低冪次m if length(x)<N x=[x,zeros(1,N-length(x))]; % 若x的長度不是2的冪，補零到2的整數冪 end nxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1; % 求1:2^m數列序號的倒序 y=x(nxd); % 將x倒序排列作為y的初始值 for mm=1:m % 將DFT作m次基2分解,從左到右，對每次分解作DFT運算,共做m級蝶形運算，每一級都有2^(mm-1)個蝶形結 Nz=2^mm;u=1; % 旋轉因子u初始化為WN^0=1 WN=exp(-i*2*pi/Nz); % 本次分解的基本DFT因子WN=exp(-i*2*pi/Nz) for j=1:Nz/2 % 本次跨越間隔內的各次蝶形運算,在進行第mm級運算時需要2^(mm-1)個 蝶形 for k=j:Nz:N % 本次蝶形運算的跨越間隔為Nz=2^mm kp=k+Nz/2; % 蝶形運算的兩個因子對應單元下標的關系 t=y(kp)*u; % 蝶形運算的乘積項 y(kp)=y(k)-t; % 蝶形運算 y(k)=y(k)+t; % 蝶形運算 end u=u*WN; % 修改旋轉因子,多乘一個基本DFT因子WN end end y y1=fft(x) %自己編的FFT跟直接調用的函數運算以後的結果進行對比
⑽ 分裂基FFT演算法和基2FFT演算法在matlab里怎麼模擬實現啊
matlab里不是直接可以用fft與2FFT的嗎？
你需要重寫嗎？那查看一下啊fft.m文件唄
%% FFT of Signals in Matlab
% August 31, 2010
% Robert Francis
% [email protected]
%
%% Define time vector
samplingFrequency = 1000; %Hz
timeStep = 1/samplingFrequency; %sec
T = 1; %sec
S = T*samplingFrequency; %samples
samples = 1:S; %samples
time = samples*timeStep; %sec
%% Generate Sine wave
frequency1 = 5; %Hz
frequency2 = 7; %Hz
sineWave = sin(time*frequency1*2*pi)+sin(time*frequency2*2*pi);
figure(1), plot(time,sineWave,'r')
xlabel('Time (sec)')
ylabel('Amplitude')
title(['Original Sinusoid with ' num2str(frequency1) 'Hz and ' num2str(frequency2) 'Hz'])
set(gca,'YLim',[-2.1 2.1])
%% Fast Fourier Transform of Sine Wave
L = length(sineWave);
NFFT = 2^nextpow2(L); % Next power of 2 from length of signal
FsineWave = fft(sineWave,NFFT)/L;
f = samplingFrequency/2*linspace(0,1,NFFT/2+1);
%% Plot single-sided amplitude spectrum.
figure, plot(f,2*abs(FsineWave(1:NFFT/2+1)))
title('Single-Sided Amplitude Spectrum of SineWave')
xlabel('Frequency (Hz)')
ylabel('|SineWave(f)|')
%% Generate square wave
pulseFrequency = 8; %Hz
squareWave = square(time*pulseFrequency*2*pi);
squareWavePos = (squareWave+1)/2;
figure, plot(time,squareWavePos,'r')
set(gca,'YLim',[-0.1 1.1])
%% Fast Fourier Transform of Square Wave
L = length(squareWave);
NFFT = 2^nextpow2(L); % Next power of 2 from length of signal
FsqWave = fft(squareWave,NFFT)/L;
f = samplingFrequency/2*linspace(0,1,NFFT/2+1);
%% Plot single-sided amplitude spectrum.
figure, plot(f,2*abs(FsqWave(1:NFFT/2+1)))
title('Single-Sided Amplitude Spectrum of SquareWave')
xlabel('Frequency (Hz)')
ylabel('|SquareWave(f)|')
%% Plot Single-sided Phase Plots
figure, plot(f,phase(FsqWave(1:NFFT/2+1)))
title('Single-Sided Phase Spectrum of SquareWave')
xlabel('Frequency (Hz)')
ylabel('Phase(Y(f))')
%% Plot Amplitude and Phase in Same figure
figure,
subplot(2,1,1), plot(f,2*abs(FsqWave(1:NFFT/2+1)))
title('Single-Sided Amplitude Spectrum of SquareWave')
xlabel('Frequency (Hz)')
ylabel('|SquareWave(f)|')
subplot(2,1,2), plot(f,phase(FsqWave(1:NFFT/2+1)))
title('Single-Sided Phase Spectrum of SquareWave')
xlabel('Frequency (Hz)')
ylabel('Phase(SquareWave(f))')
%% Decibels
magnitude = 2*abs(FsqWave(1:NFFT/2+1));
dbs = 20*log10(magnitude);
dbs_anotherWay = mag2db(magnitude);
dbs_yetAnotherWay = pow2db(magnitude); %Note: 10*log10(magnitude)
figure,
subplot(2,1,1), plot(f,dbs)
title('Single-Sided Amplitude Spectrum of SquareWave')
xlabel('Frequency (Hz)')
ylabel('|SquareWave(f)| (db)')
subplot(2,1,2), plot(f,phase(FsqWave(1:NFFT/2+1)))
title('Single-Sided Phase Spectrum of SquareWave')
xlabel('Frequency (Hz)')
ylabel('Phase(Y(f))')
%% Add noise to squareWave
noisySqW = squareWavePos + 3*randn(size(squareWavePos));
figure, plot(time,noisySqW), xlabel('Time (sec)'), ylabel('Amplitude')
%% Fourier Transform of noisySqW
L = length(noisySqW);
NFFT = 2^nextpow2(L); % Next power of 2 from length of signal
FnsqWave = fft(noisySqW,NFFT)/L;
f = samplingFrequency/2*linspace(0,1,NFFT/2+1);
%% Plot noisy Fourier Transform
figure,
subplot(2,1,1), plot(f,2*abs(FnsqWave(1:NFFT/2+1)))
title('Single-Sided Amplitude Spectrum of SquareWave')
xlabel('Frequency (Hz)')
ylabel('|SquareWave(f)|')
subplot(2,1,2), plot(f,phase(FnsqWave(1:NFFT/2+1)))
title('Single-Sided Phase Spectrum of SquareWave')
xlabel('Frequency (Hz)')
ylabel('Phase(Y(f))')
%% Fourier Transform of an Impulse
impulse = [0 0 0 0 0 1 0 0 0 0 0];
L = length(impulse);
NFFT = 2^nextpow2(L); % Next power of 2 from length of signal
Fimpulse = fft(impulse,NFFT)/L;
figure, plot(abs(Fimpulse))
%% Fourier Transform of a DC Signal
DC = ones(1,1000);
L = length(DC);
NFFT = 2^nextpow2(L); % Next power of 2 from length of signal
FDC = fft(DC,NFFT)/L;
DC2 = ones(1,100000);
L2 = length(DC2);
NFFT2 = 2^nextpow2(L2); % Next power of 2 from length of signal
FDC2 = fft(DC2,NFFT2)/L2;
% Arbitrary frequency vectors
f1 = linspace(0,1,length(FDC));
f2 = linspace(0,1,length(FDC2));
figure, plot(f1,abs(FDC),'b'), hold on, plot(f2,abs(FDC2),'r')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2FFT
%% FFT of Signals in Matlab
% August 31, 2010
% Robert Francis
% [email protected]
%
%% Define time vector
samplingFrequency = 1000; %Hz
timeStep = 1/samplingFrequency; %sec
T = 1; %sec
S = T*samplingFrequency; %samples
samples = 1:S; %samples
time = samples*timeStep; %sec
%% Generate Sine wave
frequency1 = 5; %Hz
frequency2 = 7; %Hz
sineWave = sin(time*frequency1*2*pi)+sin(time*frequency2*2*pi);
figure(1), plot(time,sineWave,'r')
xlabel('Time (sec)')
ylabel('Amplitude')
title(['Original Sinusoid with ' num2str(frequency1) 'Hz and ' num2str(frequency2) 'Hz'])
set(gca,'YLim',[-2.1 2.1])
%% Fast Fourier Transform of Sine Wave
L = length(sineWave);
NFFT = 2^nextpow2(L); % Next power of 2 from length of signal
FsineWave = fft(sineWave,NFFT)/L;
f = samplingFrequency/2*linspace(0,1,NFFT/2+1);
%% Plot single-sided amplitude spectrum.
figure, plot(f,2*abs(FsineWave(1:NFFT/2+1)))
title('Single-Sided Amplitude Spectrum of SineWave')
xlabel('Frequency (Hz)')
ylabel('|SineWave(f)|')
%% Generate square wave
pulseFrequency = 8; %Hz
squareWave = square(time*pulseFrequency*2*pi);
squareWavePos = (squareWave+1)/2;
figure, plot(time,squareWavePos,'r')
set(gca,'YLim',[-0.1 1.1])
%% Fast Fourier Transform of Square Wave
L = length(squareWave);
NFFT = 2^nextpow2(L); % Next power of 2 from length of signal
FsqWave = fft(squareWave,NFFT)/L;
f = samplingFrequency/2*linspace(0,1,NFFT/2+1);
%% Plot single-sided amplitude spectrum.
figure, plot(f,2*abs(FsqWave(1:NFFT/2+1)))
title('Single-Sided Amplitude Spectrum of SquareWave')
xlabel('Frequency (Hz)')
ylabel('|SquareWave(f)|')
%% Plot Single-sided Phase Plots
figure, plot(f,phase(FsqWave(1:NFFT/2+1)))
title('Single-Sided Phase Spectrum of SquareWave')
xlabel('Frequency (Hz)')
ylabel('Phase(Y(f))')
%% Plot Amplitude and Phase in Same figure
figure,
subplot(2,1,1), plot(f,2*abs(FsqWave(1:NFFT/2+1)))
title('Single-Sided Amplitude Spectrum of SquareWave')
xlabel('Frequency (Hz)')
ylabel('|SquareWave(f)|')
subplot(2,1,2), plot(f,phase(FsqWave(1:NFFT/2+1)))
title('Single-Sided Phase Spectrum of SquareWave')
xlabel('Frequency (Hz)')
ylabel('Phase(SquareWave(f))')
%% Decibels
magnitude = 2*abs(FsqWave(1:NFFT/2+1));
dbs = 20*log10(magnitude);
dbs_anotherWay = mag2db(magnitude);
dbs_yetAnotherWay = pow2db(magnitude); %Note: 10*log10(magnitude)
figure,
subplot(2,1,1), plot(f,dbs)
title('Single-Sided Amplitude Spectrum of SquareWave')
xlabel('Frequency (Hz)')
ylabel('|SquareWave(f)| (db)')
subplot(2,1,2), plot(f,phase(FsqWave(1:NFFT/2+1)))
title('Single-Sided Phase Spectrum of SquareWave')
xlabel('Frequency (Hz)')
ylabel('Phase(Y(f))')
%% Add noise to squareWave
noisySqW = squareWavePos + 3*randn(size(squareWavePos));
figure, plot(time,noisySqW), xlabel('Time (sec)'), ylabel('Amplitude')
%% Fourier Transform of noisySqW
L = length(noisySqW);
NFFT = 2^nextpow2(L); % Next power of 2 from length of signal
FnsqWave = fft(noisySqW,NFFT)/L;
f = samplingFrequency/2*linspace(0,1,NFFT/2+1);
%% Plot noisy Fourier Transform
figure,
subplot(2,1,1), plot(f,2*abs(FnsqWave(1:NFFT/2+1)))
title('Single-Sided Amplitude Spectrum of SquareWave')
xlabel('Frequency (Hz)')
ylabel('|SquareWave(f)|')
subplot(2,1,2), plot(f,phase(FnsqWave(1:NFFT/2+1)))
title('Single-Sided Phase Spectrum of SquareWave')
xlabel('Frequency (Hz)')
ylabel('Phase(Y(f))')
%% Fourier Transform of an Impulse
impulse = [0 0 0 0 0 1 0 0 0 0 0];
L = length(impulse);
NFFT = 2^nextpow2(L); % Next power of 2 from length of signal
Fimpulse = fft(impulse,NFFT)/L;
figure, plot(abs(Fimpulse))
%% Fourier Transform of a DC Signal
DC = ones(1,1000);
L = length(DC);
NFFT = 2^nextpow2(L); % Next power of 2 from length of signal
FDC = fft(DC,NFFT)/L;
DC2 = ones(1,100000);
L2 = length(DC2);
NFFT2 = 2^nextpow2(L2); % Next power of 2 from length of signal
FDC2 = fft(DC2,NFFT2)/L2;
% Arbitrary frequency vectors
f1 = linspace(0,1,length(FDC));
f2 = linspace(0,1,length(FDC2));
figure, plot(f1,abs(FDC),'b'), hold on, plot(f2,abs(FDC2),'r')