信号处理 #

% 参数设置
fs = 1000;              % 采样频率
t = 0:1/fs:1-1/fs;      % 时间向量

% 正弦波
f = 50;                 % 频率
sin_wave = sin(2*pi*f*t);

% 余弦波
cos_wave = cos(2*pi*f*t);

% 方波
square_wave = square(2*pi*f*t);

% 锯齿波
sawtooth_wave = sawtooth(2*pi*f*t);

% 脉冲信号
pulse = pulstran(t, [0.5], 'rectpuls', 0.1);

% 白噪声
noise = randn(size(t));

% 添加噪声
signal = sin(2*pi*50*t) + 0.5*randn(size(t));

% 线性扫频
chirp_signal = chirp(t, 0, 1, 500);  % 0到500Hz

% 对数扫频
chirp_log = chirp(t, 10, 1, 500, 'logarithmic');

% 信号
fs = 1000;
t = 0:1/fs:1-1/fs;
signal = sin(2*pi*50*t) + 0.5*sin(2*pi*120*t);

% FFT
N = length(signal);
Y = fft(signal);
P2 = abs(Y/N);
P1 = P2(1:N/2+1);
P1(2:end-1) = 2*P1(2:end-1);

% 频率轴
f = fs*(0:(N/2))/N;

% 绘图
plot(f, P1);
xlabel('频率 (Hz)');
ylabel('幅值');
title('单边频谱');

% 使用pwelch
[pxx, f] = pwelch(signal, [], [], [], fs);
plot(f, 10*log10(pxx));
xlabel('频率 (Hz)');
ylabel('功率谱密度 (dB/Hz)');

% 使用periodogram
[pxx, f] = periodogram(signal, [], [], fs);

% STFT
[s, f, t_stft] = spectrogram(signal, 256, 250, 256, fs);

% 绘制频谱图
imagesc(t_stft, f, abs(s));
axis xy;
xlabel('时间 (s)');
ylabel('频率 (Hz)');
colorbar;

% 低通滤波器
fs = 1000;
fc = 100;              % 截止频率
order = 100;           % 滤波器阶数

b = fir1(order, fc/(fs/2));  % FIR低通

% 频率响应
freqz(b, 1, 1024, fs);

% 应用滤波器
filtered = filter(b, 1, signal);

% Butterworth滤波器
order = 6;
fc = 100;

[b, a] = butter(order, fc/(fs/2));  % 低通

% 带通滤波器
f_low = 50;
f_high = 150;
[b, a] = butter(order, [f_low, f_high]/(fs/2), 'bandpass');

% 应用滤波器
filtered = filtfilt(b, a, signal);  % 零相位滤波

% Butterworth - 最平坦通带
[b, a] = butter(order, Wn);

% Chebyshev Type I - 通带波纹
[b, a] = cheby1(order, Rp, Wn);

% Chebyshev Type II - 阻带波纹
[b, a] = cheby2(order, Rs, Wn);

% Elliptic - 最陡峭过渡带
[b, a] = ellip(order, Rp, Rs, Wn);

% 设计滤波器
[b, a] = butter(6, 0.2);

% 频率响应
[h, w] = freqz(b, a, 512);

% 绘图
figure;
subplot(2, 1, 1);
plot(w/pi, 20*log10(abs(h)));
xlabel('归一化频率');
ylabel('幅值 (dB)');
title('幅频响应');

subplot(2, 1, 2);
plot(w/pi, unwrap(angle(h)));
xlabel('归一化频率');
ylabel('相位 (rad)');
title('相频响应');

% 打开滤波器设计工具
filterDesigner;

% 或使用命令行
designfilt('lowpassfir', 'FilterOrder', 100, 'CutoffFrequency', 100, ...
           'SampleRate', 1000);

% 生成含噪信号
fs = 1000;
t = 0:1/fs:1;
clean = sin(2*pi*50*t);
noisy = clean + 0.5*randn(size(t));

% 设计低通滤波器
fc = 100;
[b, a] = butter(6, fc/(fs/2));

% 滤波
filtered = filtfilt(b, a, noisy);

% 绘图
figure;
subplot(3, 1, 1);
plot(t, clean);
title('原始信号');

subplot(3, 1, 2);
plot(t, noisy);
title('含噪信号');

subplot(3, 1, 3);
plot(t, filtered);
title('滤波后信号');

% 模拟ECG信号
fs = 360;
t = 0:1/fs:10;
ecg = 0.5*sin(2*pi*1*t) + sin(2*pi*5*t) + 0.3*randn(size(t));

% 带通滤波 (0.5-40 Hz)
[b, a] = butter(4, [0.5 40]/(fs/2), 'bandpass');
ecg_filtered = filtfilt(b, a, ecg);

% 去除基线漂移
[b, a] = butter(4, 0.5/(fs/2), 'highpass');
ecg_baseline = filtfilt(b, a, ecg_filtered);