基于Matlab的信号处理:信号波形恢复与求各阶谐波及数据拟合
基于matlab的信号处理,信号波形恢复,求各阶谐波,数据拟合
最近在搞信号处理的项目,偶然发现Matlab的谐波分解功能挺有意思。今天顺手记录几个实用技巧,主要关于波形恢复、谐波提取和数据拟合。先从简单的合成信号开始吧,毕竟真实数据总是夹杂着各种噪声和干扰。
先整个带噪声的正弦波玩玩:
fs = 1000; % 采样率 t = 0:1/fs:1; % 时间轴 f0 = 50; % 基频 signal = 0.7*sin(2*pi*f0*t) + 0.3*cos(2*pi*3*f0*t); % 三阶谐波 noise = 0.2*randn(size(t)); % 高斯噪声 noisy_signal = signal + noise; % 画原始信号 figure; subplot(2,1,1) plot(t, noisy_signal) title('含噪信号波形')这时候得到的波形像被狗啃过似的毛刺。用移动平均滤波试试恢复:
window_size = 15; % 滑动窗口大小 b = (1/window_size)*ones(1, window_size); filtered_signal = filter(b, 1, noisy_signal); % 对比滤波效果 subplot(2,1,2) plot(t, filtered_signal, 'r', 'LineWidth',1.5) hold on plot(t, signal, 'k--') legend('滤波后','原始信号')这里filter函数的第一个参数是分子系数,第二个是分母系数。实际效果取决于窗口大小——太小去噪不彻底,太大会导致相位延迟。我试过窗口取采样周期的1/10左右效果比较稳。
基于matlab的信号处理,信号波形恢复,求各阶谐波,数据拟合
接下来搞谐波分解,FFT肯定是基本功:
N = length(signal); Y = fft(filtered_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; figure; stem(f,P1) xlabel('频率 (Hz)') title('单边幅度谱') axis([0 200 0 0.5])频谱图上应该能看到50Hz和150Hz两个峰。不过实际项目中会遇到频谱泄露,这时候加个汉宁窗会好很多:
window = hann(N); Y_windowed = fft(filtered_signal .* window');想单独提取各次谐波?可以玩点花样:
harmonics = cell(3,1); for k = 1:3 [~,idx] = findpeaks(P1, 'NPeaks',k,'SortStr','descend'); harmonics{k} = P1(idx(k)); end disp(['各次谐波幅度:', num2str(cell2mat(harmonics)')])不过这种方法在密集频谱中可能翻车,遇到这种情况建议改用谐波重构:
t_harmonic = 0:1/fs:0.1; % 取局部波形 recon = harmonics{1}*sin(2*pi*f0*t_harmonic) + harmonics{3}*cos(2*pi*3*f0*t_harmonic);最后说说数据拟合,polyfit虽然老套但实用:
x = t(1:50:end); % 降采样 y = real(filtered_signal(1:50:end)); p = polyfit(x,y,5); % 5次多项式 fitted_curve = polyval(p,x); figure; scatter(x,y,'filled') hold on plot(x,fitted_curve,'LineWidth',2)注意多项式阶数别太高,否则会过拟合。遇到周期性数据可以试试傅里叶级数拟合:
fit_type = fittype('fourier3'); [fitted_model, gof] = fit(x', y', fit_type);这里gof结构体里的rsquare值大于0.9基本算合格。实际项目中遇到过采样点分布不均匀导致拟合失败,这时候建议先做插值再拟合。
这些方法组合起来,基本上能应付大部分波形处理需求。不过要警惕采样定理——之前有次没注意混叠,结果重构出个四不像的波形,血泪教训啊。
