时频分析算法应用-matlab

在网上收集的点时频分析代码，可以直接运行的那种

原理大家可以自学一下，直接上代码！

1、首先输入一个时域信号，这里加了噪声，方便展示效果

clear； clc;

%% 参数设置
fs =1000;                     % 采样率
n = length(data(:,1));        % 采样点数
T = n/1000;                   % 采样时间
t = [0:1/fs:T-1/fs]';         % 采样时间向量


%% 输入一个信号，这里设计了一个雷克子波

w1.r = 10;    % 子波宽度
w1.t0 = 0.3;  % 震发时刻
w1.f =  120;  % 主频
w1.src = 2;   % 尺度
w1.wavelet11= w1.src*exp(-(2*pi*w1.f/w1.r)^2*t.^2).*sin(2*pi*w1.f.*t);
w1.wavelet = zeros(n,1);
w1.wavelet(w1.t0*fs+1:end) =  w1.wavelet11(1:end-w1.t0*fs);

wavelet = w1.wavelet;

%% 信号加噪
Noise.sigma = 0.3; 
Noise.noise1 = normrnd(0,Noise.sigma,n,1); % 随机噪声


Noise.fg1 = 50; Noise.src1 = 0.3;     % 工频干扰
Noise.fg2 = 150; Noise.src2 = 0.1;
Noise.noise2 = Noise.src1*sin(2*pi*Noise.fg1*t) + Noise.src2*sin(2*pi*Noise.fg2*t);
signal = wavelet + Noise.noise1 + Noise.noise2;
signal = signal/max(abs(signal)) ;% 归一化

SNR.signalPower = sum(abs(wavelet).^2)/length(wavelet);             % 求出信号功率
SNR.noisePower=sum(abs(signal-wavelet).^2)/length(signal);          % 求出噪声功率
SNR.SNR=SNR.signalPower/SNR.noisePower;                % 由信噪比定义求出信噪比，单位为db
disp([SNR.SNR,'   数据信噪比 :  ' , num2str(SNR.SNR)]);

2、一般当数据质量较差的情况下需要进行噪声压制，也就是信噪分离技术，常用机器学习类算法，在知道数据时频特征的前提下，也可以用先用带通滤波进行初步去噪，当然也有小波阈值去噪（但似乎没啥效果，特别是当有效信号为非稳态时），怎么简单怎么来！咱们就先带通一下，反正就两行代码的事儿

%% 带通滤波
signal_bandpass = bandpass(signal,[80 150],fs);
signal_bandpass = signal_bandpass/max(abs(signal_bandpass))  ;% 归一化

3、短时傅里叶变换

针对非稳态信号，需要分析其时频谱特性，常用的FFT不管啥用，不能知道该频段出现的时间，为此，可以用短时傅里叶变换（STFT），但受窗函数影响，精度有限，效果就那样吧

Stft.windowLength = 64; % 窗口长度
Stft.hopSize = 8;       % 帧移大小
 
[Stft.S, Stft.f, Stft.tt] = spectrogram(signal, Stft.windowLength, Stft.windowLength-Stft.hopSize, Stft.windowLength, fs, 'yaxis');

Stft.Abs_S = abs(Stft.S);

4、小波变换

有时频分析基础多多少少都知道小波变换就是从STFT发展而来的，克服了诸多缺点，处理非稳态信号效果较好，但也有时间域和频率域分辨率相互矛盾等缺陷。。。总体来说，精度还是挺不错的

%% 小波变换
WT.wavename = 'cmor1-1';
WT.totalscal = 4096;
WT.Fc = centfrq(WT.wavename); % 小波的中心频率
WT.c = 2*WT.Fc*WT.totalscal;
WT.scals = (WT.c)./(1:WT.totalscal);
f = scal2frq(WT.scals,WT.wavename,1/fs);  % 将尺度转换为频率
coefs_signal = cwt(signal,WT.scals,WT.wavename); Coefs_signal = abs(coefs_signal); Coefs_signal  = Coefs_signal/max(max(Coefs_signal));
coefs_bandpass = cwt(signal_bandpass,WT.scals,WT.wavename); Coefs_bandpass = abs(coefs_bandpass);

5、经验模态分解

这个时频分析方法可以用于信噪分离，也有很多对它改进的，如EEMD、CEEMDAN，效果大差不差。但在处理本次的非稳态信号，它们都没效果

%% 经验模态分解
[Emd.imf,Emd.residual] = emd(signal);
Emd.L_imf = length(Emd.imf(1,:));
figure(1)
for i = 1:Emd.L_imf
    subplot(Emd.L_imf,1,i)
    plot(t,Emd.imf(:,i));set(gca,'XLim',[0 2],'YLim',[-1.0 1.0],'FontSize',13)
end
signal = Emd.imf(:,1);

6、变分模态分解

咳。。。长叹息以掩涕兮，你们学吧，代码在这

%% 变分模态分解
Vmd.alpha = 2500;       %  适度的带宽约束/惩罚因子
Vmd.tau = 0;          %  噪声容限（没有严格的保真度执行）
Vmd.K = 6;              %  分解的模态数，可以自行设置，这里以8为例。
Vmd.DC = 0;             % 无直流部分
Vmd.init = 1;           %  均匀初始化
Vmd.tol = 1e-6;        
[Vmd.u, Vmd.u_hat, Vmd.omega] = VMD(signal, Vmd.alpha, Vmd.tau, Vmd.K, Vmd.DC, Vmd.init, Vmd.tol); %其中u为分解得到的IMF分量
Vmd.u = Vmd.u'; Vmd.Lu = length(Vmd.u(1,:));
figure(1)
for i = 1:Vmd.Lu
    subplot(Vmd.Lu,1,i)
    plot(t,Vmd.u(:,i))
    set(gca,'XLim',[0 2],'YLim',[-1.0 1.0],'FontSize',14 , 'FontName', 'Times')
end

7、S变换

论文里面效果一个个贼好，为啥没效果，给我个解释，我真的会谢

%% S变换
St.frelow = 0;
St.frehigh = 500;
St.alpha = 0.1;
St.fre = St.frelow:St.alpha:St.frehigh;
[St.ST,St.times,St.frequencies] = st(signal);

当然这里展示的只是一些基础的时频分析方法，还有如SST、SET等等，分析效果都有其独特的优势，这里就不一一展示了

8、成图！

最后放一个成图的代码，图好看，论文就成功了一半！下面是一个成图模板，能用

figure(1)
plot(t,wavelet,'b-','linewidth',1.5);legend('Ricker Wavelet','FontSize',24, 'FontName', 'Times');grid on
set(gca,'XLim',[0 0.5],'YLim',[-2.0 2.0],'FontSize',24, 'FontName', 'Times');set(gca,'XTick',[0:0.1:0.5],'YTick',[-2:1:2],'linewidth',1.0);
xlabel('Time(s)', 'FontSize', 24,  'FontName', 'Times'),ylabel('Normalized amplitude', 'FontSize', 24,  'FontName', 'Times')

Rectangle.recx9=0.98; Rectangle.recx10=1.03; Rectangle.recy9=-1.2; Rectangle.recy10=1.2;Rectangle.recy1=80; Rectangle.recy2=180;

figure(2)
subplot(2,1,1)
imagesc(t,f,Coefs_signal);
set(gca,'YDir','normal');set(gca,'XLim',[0 2],'YLim',[0 500],'FontSize',24, 'FontName', 'Times')
set(gca,'YTick',[0:100:500])
colorbar('FontSize',20); colormap('jet');
xlabel('Time(s)', 'FontSize', 24, 'FontName', 'Times'); ylabel('Frequency (Hz)', 'FontSize', 24, 'FontName', 'Times');
rectangle('Position',[Rectangle.recx9 Rectangle.recy1 Rectangle.recx10-Rectangle.recx9 Rectangle.recy2-Rectangle.recy1],'EdgeColor','white','LineStyle','--','Linewidth',1.5);
subplot(2,1,2)
imagesc(t,f,Coefs_bandpass);
set(gca,'YDir','normal');set(gca,'XLim',[0 2],'YLim',[0 500],'FontSize',24, 'FontName', 'Times')
set(gca,'YTick',[0:100:500])
colorbar('FontSize',20);colormap('jet');
xlabel('Time(s)', 'FontSize', 24, 'FontName', 'Times'); ylabel('Frequency (Hz)', 'FontSize', 24, 'FontName', 'Times');
rectangle('Position',[Rectangle.recx9 Rectangle.recy1 Rectangle.recx10-Rectangle.recx9 Rectangle.recy2-Rectangle.recy1],'EdgeColor','white','LineStyle','--','Linewidth',1.5);
% print(gcf,'-dpng','-r600',['原始微震信号与去噪后曲线对比']);  % 要成就成高清图！