傅里叶变化频谱图

一 . 整体示例

示例代码创建：

%%傅里叶变换频谱图
%时域分析
ts = 0:0.01:10; 
sigl = sin(2*pi*ts);%单一成分慢信号
sig2 = 5*sin(2*pi*10*ts+. 75*pi);%单一成分快信号
subplot (511) ;plot(sig1)
subplot (512) ;plot (sig2)
%多成分
sig3 = sin(2*pi*ts) +5*sin(2*pi*10*ts+.75*pi)+...
3.5*sin (2*pi*50*ts+.25*pi) ;

• 数据信号三属性：频率，幅值，相位
• 时域分析的不足，无法让我们确切地知道信号有多少成分
• 傅里叶变换可以简单理解为：无穷多个（单一成分的）正弦函数逼近任何信号
• fft幅频图绘制：

Y = fft(sign);            %对sign进行傅里叶变换
P2 = abs(Y/L);            %对纵坐标量纲即幅值还原，针对直流分量，此处abs取模长
L = length(Y)             %Y表示变换后复数向量，L为数据长度/点数
P1 = P2(1:L/2+1);         %取(L/2)+1个数据点
P1 (2:end-1) = 2*P1 (2:end-1);        %针对直流分量，此处是围绕0震荡
%对中间点能量进行还原
%因为奈奎斯特原理，去掉了一半数据，根据能量守恒需进行补充。
f = Fs*(0:(L/2))/L;                   %将纵坐标频率化，从0到Fs/2，Fs表示采样率
%因为奈奎斯特原理，只能取到采样率一半的频率，即最大频率为Fs/2(hz)
plot(f, P1)
title(' Single-Sided Amplitude Spectrum of X(t)' )
xlabel('f(Hz)')              %频率化了的横坐标
ylabel('|P1(f)|')            %此处是选取傅里叶变化的振幅幅值同还原后的频率一一对应
%傅里叶变换将时域中的信号分解成多个正弦和余弦波的频率成分。这些频率成分在频域中表示为复数形式，实部表示振幅，虚部表示相位

二 . 能量泄露

• 傅里叶变换时会对信号进行截断和周期延拓，如果不是进行截断一个周期将会造成能量泄露。
• 应使窗函数频谱的主瓣宽度应尽量窄，以获得高的频率分辨能力；旁瓣衰减应尽量大，以减少频谱拖尾。

查看不同窗函数的特性：fdatool-->FIR-->Window-->view

ts = 0:0.01:10;
sig = sin(2*pi*ts);
sig_fft = fft(sig);
L = length(sig);
sig_amp = abs(sig_fft)/L;%
sig_pl = sig_amp(1:L/2+1);
sig_pl(2:end-1) = sig_pl(2:end-1)*2;
fs = 100*(0:L/2)/L;
subplot(211)
plot(fs,sig_pl);
wind = window('hamming',L) ;
sig_win = sig'.*wind; %需要有修正系数
sig_fft = fft(sig_win);
sig_amp = abs(sig_fft)/L;%
sig_pl = sig_amp(1:L/2+1);
sig_pl(2:end-1) = sig_pl(2:end-1)*2;
fs = 100*(0:L/2)/L;
subplot(212)
plt = plot(fs,sig_win);
plt.Color(4) = .2;

三 . 时域分析指标

• struct结构体可以包含多个字段（fields），每个字段都可以存储一个值或一个数据项。每个字段都有一个名称和与之关联的值。使用struct，可以将不同类型的数据组织在一起，形成一个逻辑上相关的数据集合。每个字段可以用于存储不同类型的数据，例如数字、字符串、数组、其他结构体等。如

% 创建结构体
person.name = 'John Smith';
person.age = 30;
person.height = 1.8;

% 访问结构体字段
disp(person.name);
disp(person.age);
disp(person.height);

%创建一个最大值存放的结构体(struct)
emg_ws.Max(1) = max (stron) ;
emg_ws.Max(2) = max (wea)
emg_ws.Min(1) = min (stron) ;
emg_ws.Min(2) = min (wea) ;
%峰峰值
emg_ws. Peak2valley(1) = max(stron) - min(stron) ;
emg_ws. Peak2valley(2) = max(wea) - min (wea) ;
%均方根值，有效值
emg_ws.RMS(1) =rms(stron);
emg_ws.RMS(2) = rms(wea) :
%绝对平均值
emg_ws.ABS(1)= mean(abs(stron));
emg_ws.ABS(2)=mean(abs(wea));
%过零点
count0 = 0;
for i = 1:length(wea)-1
    if wea(i)*wea(i+1) < 0
        countO = countO + 1;
    end
end
emg_Ws.Count0(2) = count0;
%平均值
emg_ws.Mean(1) = mean(stron) ;
emg_ws.Mean(2)= mean(wea);




%%无量纲指标
%峭度，kurtosis
%{
峭度（Kurtosis）是描述随机变量概率分布形态陡峭程度的统计量。它衡量了概率分布的尾部厚度和峰值的尖锐程度，用于描述数据的非对称性和尖峰性。
正态分布的峭度为3，被称为标准峭度（excess kurtosis）。如果数据的峭度大于3，则说明数据分布比正态分布更陡峭和尖锐；如果峭度小于3，则说明数据分布比正态分布更平坦和扁平。

偏度（Skewness）是用于描述概率分布偏斜程度的统计量。它衡量了概率分布的不对称性，即数据分布在平均值周围的左右偏移程度。

偏度可以帮助我们了解数据分布的形状，特别是尾部的倾斜方向和程度。正偏（Positive skewness）表示数据分布右偏，即尾部延伸到较大的值，而负偏（Negative skewness）表示数据分布左偏，即尾部延伸到较小的值。对称分布的偏度接近于零。
%}

emg_ws.Kurtosis(1) = kurtosis(stron) ;
emg_ws.Kurtosis(2) = kurtosis(wea) ;


%偏度
emg_ws.Skewness(1) = skewness (stron);
emg_ws.Skewness(2) = skewness (wea);


%峰值因子
emg_ws.Peakfactor(1) = max(stron) /rms(stron) ;
emg_ws.Peakfactor(2) = max(wea) /rms(wea) ;

%脉冲因子
emg_ws.Pulsefactor(1) = max(stron)/mean(abs(stron)) ;
emg_ws.Pulsefactor(2) = max(wea) /mean(abs(wea)) ;

%波形指标=脉冲因子/峰值因子
emg_ws.Wave(1)= rms(stron)/mean(abs(stron)) ;
emg_ws.Wave(2)= rms(wea) /mean(abs(wea)) ;

四 . 频域分析指标

• 基于频谱图提取频域特征

%% 方法1---基于频谱图提取频域特征
%特征1，频域幅值平均值（AF_AM）
FDF.AF_AM = mean(P1);
%特征2，重心频率(AF_CF)
FDF.AF_CF = sum(f'.*P1)/sum(P1);
%特征3，均方频率（AF_MSF）
FDF.AF_MSF = sum((f.*f)'.*P1)/sum(P1);
%特征4，方差频率(AF_RMSF)
FDF.AF_RMSF = sqrt(sum((f.*f)'.*P1)/sum(P1));
%特征5，频率方差(AF_FVAR)
FDF.AF_FVAR = sum(((f-FDF.AF_CF).^2)'.*P1)/sum(P1);
%特征...

• 基于功率谱图提取频域特征

%% 方法2---基于功率谱图提取频域特征
P1S = P1.*P1;
%特征1，平均频率（PS_MNF）
FDF.PS_MNF = sum(f'.*P1S)/sum(P1S);
%特征2，中值频率（PS_MDF）
FDF.PS_MDF = [];
%特征3，总功率（PS_TP）
FDF.PS_TP = sum(P1S);
%特征4，平均功率（PS_MNP）
FDF.PS_MNP = mean(P1S);
%特征5，最大功率值所对应的频率（PS_MPF）
[mm, ii] = max(P1S);
FDF.PS_MPF = f(ii);
%特征6，低频与高频功率比值（PS_LHR）
FDF.PS_LHR = [];
%特征...

• 计算中值频率

% 假设你已经有一个信号向量 x，并且采样频率为 Fs

% 使用pwelch函数计算信号的功率谱密度估计
[pxx, f] = pwelch(x, [], [], [], Fs);

% 计算每个频率点的权重，可以选择使用功率谱pxx或其他相关的能量值
weights = pxx;

% 计算加权频率
weighted_f = f .* weights;

% 计算总权重
total_weight = sum(weights);

% 计算中值频率
median_frequency = sum(weighted_f) / total_weight;

%{
[pxx, f] = pwelch(x, [], [], [], Fs); 是使用MATLAB中的pwelch函数来计算信号的功率谱密度估计的语句。
x：表示输入的信号向量。它是需要进行功率谱估计的信号数据。
[]：表示采用默认值的参数，用于设置窗函数。
[]：表示采用默认值的参数，用于设置重叠样本数。
[]：表示采用默认值的参数，用于设置FFT长度。
Fs：表示信号的采样频率。它指定了信号的采样率，以便在计算功率谱密度估计时正确地标定频率轴。
pwelch函数是MATLAB中用于计算功率谱密度估计的函数。它使用Welch's方法，结合窗函数和重叠技术，以获得信号在频域上的能量分布情况。pwelch函数返回两个输出变量：
pxx：表示计算得到的功率谱密度估计。它是一个向量，包含了频域上的功率值。
f：表示频率向量。它是一个与功率谱密度估计pxx相对应的频率轴，以Hz为单位。
%}

五 . 采样率问题

采样率越高，滤波器越难发挥应有作用，正常衰减应达到-20db以下。，解决方法，使用降采样方法：

resample函数：

y = resample(x, p, q)
%采用多相滤波器对时间序列进行重采样，得到的序列y的长度为原来的序列x的长度的p/q倍，p和q都为正整数。此时，默认地采用使用FIR方法设计的抗混叠的低通滤波器。

y = resample(x, p, q, n)
%采用chebyshevIIR型低通滤波器对时间序列进行重采样，滤波器的长度与n成比例，n缺省值为10.

y = resample(x, p, q, n, beta)
%beta为设置低通滤波器时使用Kaiser窗的参数，缺省值为5.

y = resample(x, p, q, b)
%b为重采样过程中滤波器的系数向量。

[y, b] = resample(x, p, q)
%输出参数b为所使用的滤波器的系数向量。

decimate函数：

y = decimate(x, r)
%将时间序列x的采样频率降低为原来的1/r，即length(y)=length(x)/r。在抽取之前，默认地采用了8阶chebyshevI型低通滤波器压缩频带。

y = decimate(x, r, n)
%采用n阶chebyshevI型低通滤波器。

y = decimate(x, r, ‘fir’)
%采用30阶的FIR型低通滤波器来压缩频带，对时间序列进行整数倍抽取。

y = decimate(x, r， n, ‘fir’)
%指定当对时间序列进行整数倍抽取的时候，采用n点FIR型低通滤波器来压缩频带，对时间序列进行整数倍抽取。

downsample函数：

y = downsample(x, n)
%从第一项开始，等间隔n对x采样，得到的序列为y。

y = downsample(x, n, phase)
%从第phase+1项开始，等间隔n对x采样，得到的序列为y，而0<=phase<n.

六 .滤波器效果的仿真验证

滤波器设计完成后还可以生成Simulink模型进行仿真：按照下图中数字标号进行

1. 第一步点击左边Realize Model图标；
2. 第二步勾选“Build model using basic elements”这一项，右边四个灰色的项将自动打钩；
3. 最后点击“Realize Model”，matlab将自动生成滤波器模型，在弹出的窗口中双击模型可以观察该模型的内部结构。

 使用生成的滤波器搭建一个简单的测试模型：将两个幅值为1，频率分别为10Hz、50Hz的正弦波叠加，输入滤波器后观察滤波前后的波形。仿真时间设为1s，仿真参数中求解器类型设为固定步长，求解器选择discrete（它适用于离散无连续状态的系统），步长设为0.005s（200Hz）