[matlab]周期性信号分析
目录
信号预处理
周期性特征提取方法
频谱分析
傅里叶变换
快速傅里叶变换（FFT）
周期图法
Welch法
自相关分析
时频分析
基于模型的方法
时间序列分解
应用实例
提取信号的周期性特征是一个在信号处理领域广泛应用的技术，特别是在分析周期性信号（如机械振动、脑电信号等）时尤为重要。以下是一些常用的提取信号周期性特征的方法和步骤：
信号预处理
-
去噪：由于实际信号中常含有噪声，因此在进行周期性特征提取之前，通常需要对信号进行去噪处理，以提高信号的信噪比。去噪方法包括滤波（如低通滤波、带通滤波等）、小波变换去噪等。
-
去趋势：对于某些信号，特别是像脑电信号这样的非平稳信号，可能包含低频趋势成分。这些趋势成分会干扰周期性特征的提取，因此需要进行去趋势处理。
周期性特征提取方法
频谱分析
频谱分析是将信号从时间域转换到频率域，以观察信号的频率成分及其分布特征的方法。
傅里叶变换
傅里叶变换（Fourier Transform, FT）是一种数学工具，用于将信号从时间域（或空间域）转换到频率域。对于周期性信号，傅里叶变换将其分解为一系列正弦和余弦函数的和，这些正弦和余弦函数的频率是原始信号频率的整数倍（即谐波）。频谱图显示了这些正弦和余弦函数的振幅（或功率）随频率的变化，频谱图中的峰值对应于信号中的周期性成分。
快速傅里叶变换（FFT）
快速傅里叶变换（Fast Fourier Transform, FFT）是傅里叶变换的一种高效实现算法。FFT通过减少计算复杂度（从O(N2)降低到O(NlogN)），使得对大规模数据集的频谱分析变得可行。FFT在信号处理、图像处理、音频分析等众多领域都有广泛应用。通过FFT，我们可以快速得到信号的频谱图，进而分析信号的频率成分。
Fs = 1000; % 采样频率为1000Hz
N = length(data); % 数据长度
t = (0:N-1)/Fs; % 时间向量 % FFT分析
Y = fft(data);
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;
plot(f,P1)
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|')周期图法
周期图法是一种基于FFT的频谱估计方法，它直接通过计算信号FFT的幅值平方来估计信号的功率谱密度（Power Spectral Density, PSD）。然而，周期图法的一个主要缺点是频谱估计的方差较大，即估计结果不够平滑。这主要是因为周期图法没有考虑数据之间的统计特性，如相关性和随机性。因此，周期图法通常适用于对频谱估计精度要求不高的情况。
% 计算 Lomb-Scargle 周期图谱
[frequencies, power] = plomb(data, time);% 绘制周期图谱
figure;
plot(frequencies, power);
xlabel('Frequency');
ylabel('Power');
title('Lomb-Scargle Periodogram');
grid on;Welch法
Welch法是一种改进的频谱估计方法，旨在克服周期图法频谱估计方差大的问题。Welch法通过以下几个步骤来改进频谱估计：
- 数据分段：将原始数据分为多个重叠的段（segment）。
- 加窗：对每个段应用窗函数（如汉宁窗、海明窗等），以减少频谱泄露。
- FFT：对每个加窗后的段进行FFT，得到各段的频谱。
- 平均：将所有段的频谱进行平均，以提高频谱估计的平滑度和准确性。
Welch法通过分段加窗和平均处理，有效地降低了频谱估计的方差，提高了频谱估计的精度和平滑度。因此，Welch法在需要高精度频谱估计的场合得到了广泛应用，如通信、雷达、生物医学信号处理等领域。
% 数据 data
Fs = 1000; % 采样频率为1000Hz % Welch法频谱估计
[pxx,f] = pwelch(data,[],[],[],Fs); % 绘图
figure;
plot(f,10*log10(pxx))
title('Welch Power Spectral Density Estimate')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')自相关分析
自相关函数（Autocorrelation Function, ACF）是信号与其自身在时间轴上移动（延迟）一定量后的乘积的积分（对于连续信号）或求和（对于离散信号）。它衡量了信号在不同时间点上与其过去或未来值之间的相似度。
对于离散信号x[n]，其自相关函数Rxx[m]定义为：
其中，m是延迟量，表示信号与其自身相比的偏移量。
特性
- 对称性：自相关函数通常是偶函数，即Rxx[−m]=Rxx[m]。
- 最大值在原点：对于大多数信号，自相关函数在m=0时取得最大值，因为此时信号与其自身完全重合。
- 周期性信号的周期性：如果信号是周期性的，那么其自相关函数也将是周期性的，且周期与原信号相同。特别地，在延迟等于信号周期或其整数倍时，自相关函数会出现峰值。
% 计算自相关函数
[autocorr_values, lags] = xcorr(data, 'coeff');
% 绘制自相关函数
figure;
plot(lags, autocorr_values);
xlabel('Lag');
ylabel('Autocorrelation');
title('Autocorrelation Function');
grid on;时频分析
定义：时频分析是同时考虑信号在时间和频率域的特征，以揭示信号的时变频率特性
常用方法
- 短时傅里叶变换（STFT）：将信号划分为多个短时窗，对每个窗内的信号进行傅里叶变换，从而得到信号随时间变化的频谱。
- 小波变换：通过选择合适的小波基函数，对信号进行多尺度分析，以揭示信号在不同时间尺度上的频率特性。
- 希尔伯特-黄变换（HHT）：包括经验模态分解（EMD）和希尔伯特谱分析，适用于非线性、非平稳信号的分析
y=data
imf=emd(y)%进行EMD分解各个固有模态函数IMF（i）% 绘制原始信号
figure;
plot(t, y); % 绘制原始信号
title('原始信号');
xlabel('时间 (s)'); %时间单位秒
ylabel('幅值');
grid on; figure;
for i = 1:size(imf, 2) % 遍历 imf 的列subplot(size(imf, 2)+1, 1, i+1); plot(t, imf(:,i)); % 绘制第 i 个 IMF title(['IMF ', num2str(i)]); xlabel('时间 (s)'); ylabel('幅值');
end基于模型的方法
定义：通过建立信号的数学模型来提取周期性特征。
常用模型：
- 自回归模型（AR）：将信号表示为自身过去值的线性组合加上噪声。
- 滑动平均模型（MA）：将信号表示为噪声的线性组合，其中噪声是过去某个时刻的噪声值。
- 自回归滑动平均模型（ARMA）：结合了AR和MA的特点，是更一般的信号模型。
时间序列分解
时间序列分解是一种将时间序列数据分解成趋势、季节性和随机成分的方法。这种方法通过去除趋势和季节性成分，使得周期性特征更加突出。具体步骤如下：
- 去除趋势：首先，将数据中的趋势成分去除，得到去趋势序列。这有助于更清晰地观察数据的周期性变化。
- 去除季节性：接着，将去趋势序列中的季节性成分去除，得到季节性序列。这一步进一步剥离了影响周期性的其他因素。
- 分析随机成分：最后，将季节性序列中的随机成分去除，得到随机序列。虽然这一步不直接用于周期性分析，但它有助于理解数据中不可预测的波动。
时间序列分解的结果可以用于预测未来的发展趋势，通常结合回归分析、指数平滑和ARIMA模型等方法进行。
% 定义周期（例如，季节性周期为12）
period = 12;% 使用 MATLAB 的时序分解函数进行季节性分解
decomp = seasdecomp(data, period);% 绘制分解结果
figure;
subplot(4, 1, 1);
plot(t, data);
title('Original Data');
xlabel('Time');
ylabel('Value');subplot(4, 1, 2);
plot(t, decomp.Trend);
title('Trend');
xlabel('Time');
ylabel('Value');subplot(4, 1, 3);
plot(t, decomp.Season);
title('Seasonal');
xlabel('Time');
ylabel('Value');subplot(4, 1, 4);
plot(t, decomp.Residual);
title('Residual');
xlabel('Time');
ylabel('Value');应用实例
- 机械振动信号：通过频谱分析和自相关分析可以提取轴承、齿轮等旋转机械部件的故障特征，如周期性脉冲冲击信号。
- 脑电信号：通过频谱分析和时频分析可以提取脑电信号的节律性特征，如α波、β波等，这些特征对于脑功能研究和疾病诊断具有重要意义。