工程中对时域信号进行处理
0 整个系统的流程框图
1 中值滤波
y=medfilt1(data,300); %中值滤波,滤除高频噪声
subplot(312);
plot(t,y);title('中值滤波后');
2 基线校准
%% --------------多项式拟合--基线校准------------------------------
opol = 6;
[p,s,mu] = polyfit(t,y,opol);
f_y = polyval(p,t,[],mu);
dt_y = y - f_y;
%将基线在0处的矫正后的数据迁移回原始数据附近(最终的目的相当于只在原始数据附近进行了基线矫正)
ZX = mean(f_y); %获得多项式值序列的均值
%基线在0处的,基线矫正后的数据与原始数据对应时间点,在纵坐标上差异(迁移)的均值
dt_data = dt_y+ZX; %将基线在0处的矫正后的数据整体迁移回原始数据附近
%在同一个figure中画出原始数据以及矫正后的数据的图像,便于观察矫正效果
subplot(313)
% plot(t,y,'r*'); %用红色※号线画出原始数据
% hold on;
plot(t,dt_data,'r-'); %用蓝色线段画出基线在原始数据附近的,矫正后的数据
hold off;
title('基线校准后');
3.1 分区间峰值扫描
%% -------------------------查找波峰点-----------------------------
N=length(dt_data);
M=150; % 设置查找的分段数
dlength=floor(N/M); % 将螺纹轮廓波形等分成 M 段, 每一段的长度
for k=1:M
DATA{k}=data((k-1)*dlength+1:k*dlength); % 保存每一段的 x 值
end
for k=1:M
[bofeng(k),index1]=max(DATA{k}); % 寻找每一段里 y 的最大值
bofeng(k)=DATA{k}(index1); % 得到与波峰对应的横坐标 x
end
figure
plot(bofeng);
xlabel('划分区域的个数');
ylabel('波峰峰值');
title('在固定区间内查找峰值点');
4.1 峰值统计
%% ------门限判定,统计在所有划分区域内大于门限的所有波峰数目-----------
bofeng_door11=bofeng>2.5;
bofeng_door12=bofeng>4;
bofeng_door1=sum(sum(bofeng_door11)); %大于2.5的峰值,进行统计,认为是正常的脉搏波
bofeng_door2=sum(sum(bofeng_door12)); %大于4的峰值,进行统计,认为是异常的脉搏波
bofengnum=[bofeng_door11;bofeng_door12];
figure
bar(bofengnum,'stacked');
3.2&4.2 小波分解&小波重构
%% ---------------------小波分解与重构---------------------
[c,l]=wavedec(dt_data,3,'db2'); %小波分解
approx=appcoef(c,l,'db2');
[cd1,cd2,cd3]=detcoef(c,l,[1,2,3]); %细节系数
figure
subplot(2,2,1)
plot(approx)
title('Approximation Coefficients-近似系数')
subplot(2,2,2)
plot(cd3)
title('Level 3 Detail Coefficients-3级细节系数')
subplot(2,2,3)
plot(cd2)
title('Level 2 Detail Coefficients-2级细节系数')
subplot(2,2,4)
plot(cd1)
title('Level 1 Detail Coefficients-1级细节系数')
x=waverec(c,l,'db2'); %小波重建函数
figure
plot(x)
title('小波重建函数');
5 时频分析
本应是用小波进行时域的功率谱分析的,但是博主真的忘了怎么弄了,只能用STFT代替下,毕竟看的是功率。
%-------------STFT时频谱-------------
[S,F,T,P1]=spectrogram(x,window1,noverlap,nfft,Fs);
P1=10*log10(P1);
figure
% subplot(5,1,1)
surf(T,F,P1,'edgecolor','none'); axis tight;
view(0,90);
xlabel('Time (Seconds)'); ylabel('Hz');title('Time-frequency Analyse');
6 归一化功率谱密度
%------------------------绘制归一化功率谱---------------------
nMax = max(plot_Pxx1); nMin = min(plot_Pxx1);
ax_norm = (plot_Pxx1 - nMin)/(nMax - nMin);
figure
plot(ax_norm);
xlabel('频率(Hz)')
ylabel('归一化功率谱密度(dB)')
title('归一化小波功率谱')
grid on
本次的程序介绍就到这里了,写的不好的地方,欢迎指正,欢迎私信或评论讨论。