[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]定义为：

$R_{xx}[m]=\sum\limits_{n=-\infty}^\infty x[n]\cdot x[n+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');