基於matlab的傅里葉變換

原文出處https://blog.csdn.net/qq_37366291/article/details/79832886

例子1
作用：使用傅里葉變換找出隱藏在噪聲中的信號的頻率成分。（指定信號的參數，采樣頻率為1 kHz，信號持續時間為1秒。）

Fs = 1000;            % 采樣頻率
T = 1/Fs;             % 采樣周期
L = 1000;             % 信號長度
t = (0:L-1)*T;        % 時間向量

%%形成一個信號，包含振幅為0.7的50hz正弦信號和振幅為1的120hz正弦信號。
S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
X = S + 2*randn(size(t));             %用零均值的白噪聲破壞信號，方差為4。


plot(1000*t(1:50),X(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('t (milliseconds)')
ylabel('X(t)')1234567891011121314

由上圖可知：從時域中我們很難觀察到信號的頻率成分。怎么辦呢？當然使用強大的傅里葉變換。

Y = fft(X);     %計算傅里葉變換,X是加噪后的信號

%%
%計算雙邊譜P2。然后計算基於P2的單面譜P1和偶值信號長度L。(不太理解。。。)
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

%%
%定義頻率域f並繪制單面振幅譜P1。由於增加的噪音，振幅不完全是0.7和1。平均而言，較長的信號產生更好的頻率近似。
f = Fs*(0:(L/2))/L;
plot(f,P1)
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|')123456789101112131415

%%
%現在，對原始的，未被損壞的信號進行傅里葉變換，並得到准確的振幅，0.7和1.0。
Y = fft(S);   %S時原始的，沒有加噪的信號。
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

plot(f,P1) 
title('Single-Sided Amplitude Spectrum of S(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|')1234567891011

加上一點自己的理解。

例子2
作用：利用傅里葉變換，將高斯脈沖從時域轉換為頻域。

Fs = 100;           % Sampling frequency
t = -0.5:1/Fs:0.5;  % Time vector
L = length(t);      % Signal length

X = 1/(4*sqrt(2*pi*0.01))*(exp(-t.^2/(2*0.01)));

plot(t,X)
title('Gaussian Pulse in Time Domain')
xlabel('Time (t)')
ylabel('X(t)')12345678910

%%
%要使用fft函數將信號轉換為頻域，首先要確定一個新的輸入長度，該輸入長度是原信號長度的下一個2次方。
%為了提高fft的性能，這將使信號X以尾隨零的形式出現。

n = 2^nextpow2(L);
Y = fft(X,n);
f = Fs*(0:(n/2))/n;
P = abs(Y/n);

plot(f,P(1:n/2+1))
title('Gaussian Pulse in Frequency Domain')
xlabel('Frequency (f)')
ylabel('|P(f)|')12345678910111213

例子3余弦波
比較時域和頻域的余弦波。指定信號的參數，采樣頻率為1kHz，信號持續時間為1秒。

Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sampling period
L = 1000;                     % Length of signal
t = (0:L-1)*T;                % Time vector

x1 = cos(2*pi*50*t);          % First row wave
x2 = cos(2*pi*150*t);         % Second row wave
x3 = cos(2*pi*300*t);         % Third row wave

X = [x1; x2; x3];

for i = 1:3
    subplot(3,1,i)
    plot(t(1:100),X(i,1:100))
    title(['Row ',num2str(i),' in the Time Domain'])
end12345678910111213141516

n = 2^nextpow2(L);
dim = 2;
Y = fft(X,n,dim);
P2 = abs(Y/n);
P1 = P2(:,1:n/2+1);
P1(:,2:end-1) = 2*P1(:,2:end-1);
for i=1:3
    subplot(3,1,i)
    plot(0:(Fs/n):(Fs/2-Fs/n),P1(i,1:n/2))
    title(['Row ',num2str(i), ' in the Frequency Domain'])
end1234567891011