为了让大学生活充实一点,多学点东西,我选修了《数字信号处理》。现在充实得不要不要的。
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)|')