经验贴：MATLAB离散序列数据处理常用手段：三次样条插值（重采样），基于FFT提取特定频率的幅值相位与信号重构，包络线实现与作用。附带：柱状图与表格绘制

三次样条插值（重采样）

插值的目的是以新的采样频率对原信号进行重构，根据不同的波形和应用场景可以选择不同的插值方法以最优还原原始信号。这里以温度信号T的三次样条插值为例，三次样条插值就是把原始数据（序列长度N）分为很N-1小区间，每个区间以下面这个公式组成一段分段函数，这样就把离散点串联起来了。通过在区间内带入自变量，即可通过插值函数得到对应因变量，从而以插值函数的形式扩展或者缩短序列。

其中
工程应用：某泵故障诊断。基于仿真模型的先验算法开发。由于算法将应用于单片机，故算法需尽量简单，数据时间分辨率能达到目的即可。首先将原始高分辨率数据重采样到每个循环140个点。

size_CSV=size(csvread('XXX.csv'),1) ; % XXX为csv文件格式，第一列为时间（角度），第二列为压力数据。
t=csvread('XXX.csv',0,0,[0,0,size_CSV-1,0]); % 从0开始计数，0,0,代表读1,1右下，后者代表另一点左上（包含该点） 
p=csvread('XXX.csv',0,1,[0,1,size_CSV-1,1]);
t=t-t(1);
spc=140; % 每转采样频率
speed=250; % 转速
tt=t0:360/spc:t0+(n+1)*360; % 对角度序列进行360/spc插值,进行低分辨率重采样
pp=interp1(t,p,tt,'spline');
p=pp; % 插值后压力
t=tt; % 插值后转角

基于FFT提取特定频率的幅值相位与信号重构

(以下代码结合了工程运用，所以部分繁琐注解希望谅解)

%需要可以拷走原始序列
p=[342.461626250000
341.456016638826
339.865479030488
337.659974486075
334.813740391819
331.305388535345
327.128733981791
328.048808750000
330.456642916273
332.777629611442
334.944160870539
336.891083738643
338.555748815510
339.878421217657
340.801761250000
341.271573016903
341.236812874543
340.649777836481
339.466354178962
337.646200724307
335.126641118416
336.614927500000
338.517347034975
340.158897190470
341.481775579691
342.431880166303
342.958879964918
343.016314065909
342.561663750000
341.556472066875
339.966406003216
337.761378656892
334.915562896668
331.407589899759
327.230843268041
328.132746250000
330.537883600631
332.856043528678
335.019611479390
336.963433243725
338.624859632209
339.944151971832
340.863980000000
341.330156185924
341.291639572688
340.700755036082
339.513388583810
337.689229108842
335.166211590552
336.645945000000
338.543604277668
340.180357118163
341.498419844408
342.443717112407
342.965943042122
343.018653958868
342.559385000000
341.549656305662
339.955053526323
337.745537390638
334.895288778509
331.382984547841
327.202457761270
328.118295000000
330.522358985813
332.839555734206
335.002254273849
336.945345174536
338.606156552537
339.924950845698
340.844390000000
341.310293815431
341.271616245802
340.680673559675
339.493356125945
337.669334782128
335.145468691836
336.633013750000
338.532118448061
340.170407392782
341.490090339993
342.437087626672
342.961073383472
343.015603687618
342.558182500000
341.550359981420
339.957818932772
337.750472542158
334.902544879036
331.392627542459
327.214408608389
328.128427500000
330.534376734521
332.853512659859
335.018219885201
336.963334195420
338.626195937375
339.947057431415
340.868577500000
341.336549865477
341.299926767888
340.711004489070
339.525655940845
337.703544843482
335.183669664765
336.657137500000
338.555651968574
340.193217326585
341.512041672206
342.458037938705
342.980890864365
343.034144509656
342.575303750000
341.565919804910
339.971685546851
337.762522863462
334.912658781203
331.400699806268
327.220177754713
328.124586250000
330.527067884854
332.842675809391
335.003804218607
336.945310959638
338.604549860013
339.921793590470
340.839718750000
341.304148485359
341.264058265388
340.671778109620
339.483209507177
337.658047854280
335.133564455287
336.617532500000
338.515079052402
340.151845580373
341.470035395511
342.415580626985
342.938167411857
342.991345870539
342.532655000000];

% p1为故障压力，感兴趣的朋友可以在后续算法自行替换
p1=[340.901917500000
339.972070353211
338.454178853652
336.318360478889
333.539009959190
330.094875085586
325.957136093356
324.529201250000
323.659951732021
322.864951287319
322.068984970311
321.200599476158
320.192040432360
318.979547373460
317.502742500000
315.705288043356
313.534779996907
310.942867992973
307.885381651907
304.322478787125
300.192986629268
300.104906250000
300.477348525783
300.662653873100
300.603179322757
300.245594022700
299.540696955032
298.443483192063
296.913937500000
294.915995875885
292.417970312470
289.392577442313
285.816986763813
281.672803467788
276.956967348088
277.256967500000
278.996377265086
280.664468508415
282.196012834600
283.528770823386
284.603457875816
285.364155905669
285.757670000000
285.734265610612
285.247629030142
284.254983789804
282.717309324916
280.599494761350
277.845529980043
279.011538750000
280.552622672218
281.842818166807
282.828790955432
283.461186770385
283.694508336756
283.487251146775
282.802071250000
281.605690324776
279.868989229638
277.567050169282
274.679200314743
271.189062507610
267.096320409319
267.993438750000
270.251677672333
272.408768048755
274.401629083485
276.169986348439
277.656376396359
278.806535417880
279.568797500000
279.894831825417
279.731744747603
283.751765282934
288.149891448785
292.447206617772
296.579107323715
300.483910000000
304.102931180991
307.380522313417
310.264182862906
312.704638991906
314.655968595199
316.075723283939
316.925062500000
317.168864230844
316.775853585212
315.718712850863
313.974206403851
311.523249222264
308.362593243545
310.231301250000
313.487808715053
316.615735878916
319.550303031162
322.228969529971
324.591507704588
326.580219888985
328.139992500000
329.218384363429
329.766033836428
329.736785385757
329.087921390900
327.780410323478
325.754197352745
327.681875000000
330.004711743264
332.047005881331
333.752224438392
335.067412890165
335.943310297774
336.334414637300
336.199111250000
335.499749163378
334.202750308892
332.278704125716
329.702455008900
326.453168513190
322.525104958112
323.665698750000
326.279431020033
328.795974724583
331.148447187277
333.272319664007
335.105547301251
336.588942763519
337.665676250000
338.282020309463
338.387367385029
337.934427506436
336.879455994734
335.182458932082
332.780951127244
334.378038750000
336.381315848195
338.118939996988
339.533437433851
340.571009147218
341.181614364375
341.319028756071
340.941023750000];

speed=250         % 曲轴转速
spc=140           % spc=sample/cycle
n1=0:360/spc:360; % 角度域横坐标
Fs=speed/60*spc;  % 信号采样频率，注意FFT进行的基础是时域信号而非角度域信号，角度域信号（n1）只为后续画图更好理解
Ts=1/Fs;          % 采样时间间隔
N=141;            % 采样信号的长度
t=(0:1:N-1)*Ts;   % 定义信号采样的时间点 t
t=t';             % 为了方便查看, 将行向量 t 转置成列向量
x=p; 			  % 定义时域采样信号 x，感兴趣的朋友可以在此替换故障压力数据

% 对时域采样信号, 执行快速傅里叶变换 FFT
X=fft(x);           % 进行离散傅里叶计算, 结果保存在 X 里
X(abs(X)<1e-8)=0;   % 将频域序列 X 中, 幅值小于 1e-8 的数值置零

% 修正频域序列的幅值, 使得 FFT 变换的结果有明确的物理意义
X=X/N;              % 将频域序列 X 除以序列的长度 N

% 提取 X 里正频率的部分, 并且将 X 里负频率的部分合并到正频率  
Y=X(1:N/2+1);              % 提取 X 里正频率的部分
Y(2:end-1)=2*Y(2:end-1);   % 将 X 里负频率的部分合并到正频率

% 计算频域序列 Y 的幅值和相角
A=abs(Y);           % 计算频域序列 Y 的幅值
Pha=angle(Y);       % 计算频域序列 Y 的相角 (弧度制)
% P=A1*cos(2*Pi*f1*(x+Pha1))+A2*cos(2*Pi*f2*(x+Pha2))+...An*cos(2*Pi*fn*(x+Phan))

% 定义序列 Y 对应的频率刻度
df=Fs/N;            % 频率间隔
f=(0:1:N/2)*df;     % 频率刻度
f=f';               % 为了方便查看, 将行向量 f 转置成列向量

% 绘制频域序列 Y 的幅频图 & 相频图 
figure(1)

subplot(3,2,[1,2]);
plot(n1,p);                    % 绘制初始角度域信号，有兴趣的朋友可以替换p1
xlabel('角度/deg.')
ylabel('初始压力/bar')

subplot(3,2,3)
stem(f(2:11),A(2:11))          % 绘制1-10倍频幅频图
xlabel('频率/Hz')
ylabel('幅值/bar')

subplot(3,2,4)
Pha=Pha/3.14159265*180 		  % 将弧度变为角度
stem(f(2:11),Pha(2:11))       % 绘制1-10倍频相频图
xlabel('频率/Hz')
ylabel('相角/deg.')

% 绘特定频率（此处仅重构1倍频）重构信号（角度域）
X=X*N       %为了IFFT准确还原时域信号
X(3:end)=0; %去除除了一倍频以外所有频率
X(1)=0;     %去除除了一倍频以外所有频率
IX = ifft(X);

subplot(3,2,[5,6])
plot(n1,IX+mean(p))  %有兴趣朋友可以替换p1
xlabel('转角/deg.')
ylabel('重构压力/bar')

下图分别是正常工况和故障工况下的频域分析展示，可以发现故障工况一倍频明显上升，且1倍频重构信号明显幅值上升（故障数据即在算法中将p替换p1，共有三处）

包络线实现与作用

包络线的主要作用是关注信号平台的变化，关注整体波动趋势而过滤小幅度高频信号的波动。在信号处理中的类似于低通滤波。此处仅实现上包络线，并且仅用了最大压力的连线。

代码如下，一定要注意序列是行还是列展开的，要不然size函数很容易出错（同样使用上一节的p1序列当纵坐标和n1序列当横坐标）

c=zeros(size(p1));
d=zeros(size(p1));
s=zeros(size(p1));
decrease=zeros(size(p1)); 
increase=zeros(size(p1)); 
speed=250         % 曲轴转速
spc=140           % spc=sample/cycle
n1=0:360/spc:360
n1=n1'

dp=diff(p1,1)/(360/spc) % 差分，360/spc为根据横坐标间隔放缩差分值（spc=140）
for i=1:size(dp,1)-1
	if dp(i)*dp(i+1)<0
		c(i)=n1(i+1); %寻找驻点横坐标
		d(i)=p1(i+1);  %寻找驻点纵坐标
	end
end 
c(find(c==0))=[]; 
d(find(d==0))=[];

figure(2)
plot(n1,p1); % 绘制原始数据
hold on;

for i=1:size(c,1)
	s(i)=find(n1==c(i)); %驻点横坐标（角度）的序列号 
end
s(find(s==0))=[]; % 序列更紧凑，下同
j=1
for i=1:size(c,1)-1
	if p1(s(i+1))>p1(s(i)) %寻找压力下降区间（s,decrease为序列号） 
    	decrease(j)=s(i); %i~i-1为下降区间（s,decrease为序列号）
        increase(j)=s(i+1); %i+1~i为上升区间（s,increase为序列号）
		j=j+1 %i~i+1为上升区间（序列号） 
	end
end
decrease(find(decrease==0))=[]; 
increase(find(increase==0))=[]; 

peak_n1=n1(increase);
peak_p1=p1(increase)
% 由于寻找极大值会漏掉原始数据的第一个和最后一个点，从而造成包络线在初末位置不完整，所以需要进行额外处理。此处由于p1是周期序列，所以对原始序列p1最后一个值赋值到peak_p1序列的第一个值和最后一个值即可。
peak_p1(2:size(peak_p1,1)+1)=peak_p1(1:end)
peak_p1(1)=p1(end)
peak_p1(size(peak_p1,1)+1)=p1(end)
% n1 同理
peak_n1(2:size(peak_n1,1)+1)=peak_n1(1:end)
peak_n1(1)=0
peak_n1(size(peak_n1,1)+1)=n1(end)
scatter(peak_n1,peak_p1) % 绘制极大值点
%对新的数据进行插值，此处三次样条将极大值序列插值回原始数据长度
peak_pp=interp1(peak_n1,peak_p1,n1,'spline');
peak_p1=peak_pp;
plot(n1,peak_p1)

柱状图和表格绘制

同样使用p1序列，直接贴代码

%% 各缸吸排液冲程平均压力直方图
% 180度判断均值 
delay=70; % 70个点可以cover180°，70=180/（360/spc）
% 后续求均值用到 sucmax_degree=[81,153,225,297,9]dramax_degree=[279,351,63,135,207]
% 根据采样率进行缩放在spc=140时，序号进行360/140的缩放
sucmax_degree=[1,28,56,84,112];%定义各缸吸液速度最大时对应角度的序号，活塞运动速度非标准正弦
dramax_degree=[70,98,126,14,42];%定义各缸吸液速度最大时对应角度的序号，活塞运动速度非标准正弦
p_all_avg=p1
Avg_P_suc1_delay=(sum(p_all_avg(sucmax_degree(1):sucmax_degree(1)+delay))/(delay+1));%计算1缸最大吸液速度时一定滞后角度内压力均值 
Avg_P_suc2_delay=(sum(p_all_avg(sucmax_degree(2):sucmax_degree(2)+delay))/(delay+1));%计算5缸最大吸液速度时一定滞后角度内压力均值
Avg_P_suc3_delay=(sum(p_all_avg(sucmax_degree(3):sucmax_degree(3)+delay))/(delay+1));%计算2缸最大吸液速度时一定滞后角度内压力均值
Avg_P_suc4_delay=((sum(p_all_avg(sucmax_degree(4):141))+sum(p_all_avg(1:(delay-(141-sucmax_degree(4))))))/(delay+1));
Avg_P_suc5_delay=((sum(p_all_avg(sucmax_degree(5):141))+sum(p_all_avg(1:(delay-(141-sucmax_degree(5))))))/(delay+1));
Avg_P_dra1_delay=(sum(p_all_avg(dramax_degree(1):dramax_degree(1)+delay))/(delay+1));%计算1缸最大吸液速度时一定滞后角度内压力均值 
Avg_P_dra2_delay=((sum(p_all_avg(dramax_degree(2):141))+sum(p_all_avg(1:(delay-(141-dramax_degree(2))))))/(delay+1));%计算5缸最大吸液速度时一定滞后角度内压力均值
Avg_P_dra3_delay=((sum(p_all_avg(dramax_degree(3):141))+sum(p_all_avg(1:(delay-(141-dramax_degree(3))))))/(delay+1));
Avg_P_dra4_delay=(sum(p_all_avg(dramax_degree(4):dramax_degree(4)+delay))/(delay+1));%计算3缸最大吸液速度时一定滞后角度内压力均值
Avg_P_dra5_delay=(sum(p_all_avg(dramax_degree(5):dramax_degree(5)+delay))/(delay+1));%计算4缸最大吸液速度时一定滞后角度内压力均值
Avg_P_delay=[Avg_P_suc1_delay,Avg_P_suc2_delay,Avg_P_suc3_delay,Avg_P_suc4_delay,Avg_P_suc5_delay,Avg_P_dra1_delay,Avg_P_dra2_delay,Avg_P_dra3_delay,Avg_P_dra4_delay,Avg_P_dra5_delay]; 


subplot(1,2,1); % 绘制吸引平均压力直方图
title(['循环各缸吸液平均压力'],'Fontsize',20);
hold on;
bar(Avg_P_delay(1:5));
set(gca,'XTickLabel',{'1','2','3','4','5'});
for i=1:5
text(i,Avg_P_delay(i)-0.03,num2str(Avg_P_delay(i)),'VerticalAlignment','bottom','HorizontalAlignment','center');%就是用test加数值，这个0.03看情况定，根据数值大小，再改就好了
end
ylim auto;
ylabel('Avg_P:Suction');

subplot(1,2,2); % 绘制排液平均压力直方图
title(['循环各缸排液平均压力'],'Fontsize',20);
hold on;
bar(Avg_P_delay(6:10));
set(gca,'XTickLabel',{'1','2','3','4','5'});
for i=6:10 
% 要注意test第一项为i-5，要不然数据对不齐
text(i-5,Avg_P_delay(i)-0.03,num2str(Avg_P_delay(i)),'VerticalAlignment','bottom','HorizontalAlignment','center');%就是用test加数值，这个0.03看情况定，根据数值大小，再改就好了
end
ylim auto;
ylabel('Avg_P:Drainage');

%% 1-10倍频表格
% 先求得频率（同上FFT）
speed=250         % 曲轴转速
spc=140           % spc=sample/cycle
n1=0:360/spc:360; % 角度域横坐标
Fs=speed/60*spc;  % 信号采样频率，注意FFT进行的基础是时域信号而非角度域信号，角度域信号（n1）只为后续画图更好理解
Ts=1/Fs;          % 采样时间间隔
N=141;            % 采样信号的长度
t=(0:1:N-1)*Ts;   % 定义信号采样的时间点 t
t=t';             % 为了方便查看, 将行向量 t 转置成列向量
x=p1; 			  % 定义时域采样信号 x，感兴趣的朋友可以在此替换故障压力数据

% 对时域采样信号, 执行快速傅里叶变换 FFT
X=fft(x);           % 进行离散傅里叶计算, 结果保存在 X 里
X(abs(X)<1e-8)=0;   % 将频域序列 X 中, 幅值小于 1e-8 的数值置零

% 修正频域序列的幅值, 使得 FFT 变换的结果有明确的物理意义
X=X/N;              % 将频域序列 X 除以序列的长度 N

% 提取 X 里正频率的部分, 并且将 X 里负频率的部分合并到正频率  
Y=X(1:N/2+1);              % 提取 X 里正频率的部分
Y(2:end-1)=2*Y(2:end-1);   % 将 X 里负频率的部分合并到正频率

% 计算频域序列 Y 的幅值和相角
A=abs(Y);           % 计算频域序列 Y 的幅值
Pha=angle(Y);       % 计算频域序列 Y 的相角 (弧度制)
% P=A1*cos(2*Pi*f1*(x+Pha1))+A2*cos(2*Pi*f2*(x+Pha2))+...An*cos(2*Pi*fn*(x+Phan))

% 定义序列 Y 对应的频率刻度
df=Fs/N;            % 频率间隔
f=(0:1:N/2)*df;     % 频率刻度
f=f';               % 为了方便查看, 将行向量 f 转置成列向量

data_row1=A(2:11);
data_row2=Pha(2:11);
data=[data_row1,data_row2];

%% 生成表格行列名称，m行n列
str1='幅值';str2='相位';str0='倍频'
n=10;
row_name=strcat(num2str((1:n)'),str0);
column_name={str1,str2};

%% 表格作图
set(figure(10),'position',[200 200 450 330]);
uitable(gcf,'Data',data,'Position',[20 20 400 290],'Columnname',column_name,'Rowname',row_name);

n1=n1'
p_all_avg=p_all_avg'

------------感谢浏览到这里------------