為了讓大學生活充實一點，多學點東西，我選修了《數字信號處理》。現在充實得不要不要的。

clc
close all
clear%=========參數設置=========%
Fs = 1000;            % Sampling frequency                    
T = 1/Fs;               % Sampling period       
L = 1500;              % Length of signal
t = (0:L-1)*T;        % Time vector
S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);     % original signal
X = S + 2*randn(size(t));     % Corrupt the signal with zero-mean white noise with a variance of 4%======畫出加噪聲的信號======%
figure
plot(1000*t(1:50),X(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('t (milliseconds)')
ylabel('X(t)')%======原始信號雙邊帶頻譜=======%
amplitude = abs(fftshift(fft(S, L)))/L;    % 計算幅度譜，fftshift是為了交換正負頻率譜，fft后緊接著要除信號長度
f = (0:L-1)/L*Fs - Fs/2 ; % 計算頻率軸，減Fs/2是為了正確顯示0頻位置
figure
plot(f, amplitude)
xlabel('頻率(Hz)')
ylabel('幅度')
title('原始信號雙邊帶頻譜')%======帶有噪聲的信號雙邊帶頻譜=======%
amplitude = abs(fftshift(fft(X, L)))/L;
f = (0:L-1)/L*Fs - Fs/2;
figure
plot(f, amplitude)
xlabel('頻率(Hz)')
ylabel('幅度')
title('帶有噪聲的信號雙邊帶頻譜')%======原始信號雙邊帶頻譜(補零)=======%
L1 = 2^nextpow2(L); % 計算FFT長度
amplitude = abs(fftshift(fft(S, L1)))/L;    % 計算幅度譜，fftshift是為了交換正負頻率譜，fft后緊接著要除信號長度
f = (0:L1-1)/L1*Fs - Fs/2 ; % 計算頻率軸，長度要用fft的長度而不是信號長度，減Fs/2是為了正確顯示0頻位置
figure
plot(f, amplitude)
xlabel('頻率(Hz)')
ylabel('幅度')
title('原始信號雙邊帶頻譜(補零)')%======帶有噪聲的信號雙邊帶頻譜(補零)=======%
L1 = 2^nextpow2(L);
amplitude = abs(fftshift(fft(X, L1)))/L;
f = (0:L1-1)/L1*Fs - Fs/2;  % 計算頻率軸，最大原信號最大頻率不超過采樣頻率的一半
figure
plot(f, amplitude)
xlabel('頻率(Hz)')
ylabel('幅度')
title('帶有噪聲的信號雙邊帶頻譜(補零)')%======加噪聲信號單邊頻譜======% 
Y = fft(X);
P2 = abs(Y/L);  % 計算幅度譜
P1 = P2(1:L/2+1);   % 提取正頻率部分，不知道為啥多取一個點
P1(2:end-1) = 2*P1(2:end-1);    % 實信號頻譜對稱
f = Fs*(0:(L/2))/L; 
figure
plot(f,P1) 
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|')%======原始信號單邊頻譜======%
Y = fft(S);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
figure
plot(f,P1) 
title('Single-Sided Amplitude Spectrum of S(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|')