小波相干性显著性检验（MATLAB R2018A）

交叉小波常被用于检测不同信号之间的相关性，其在时频域建立了不同信号之间的联系。对于两个时域信号，其交叉小波变换和交叉小波尺度谱如下：

以轴承振动信号为例，利用正常轴承与故障轴承的振动信号、故障轴承和故障轴承的振动信号分别作连续交叉小波分析得到正常-故障信号小波相干谱、故障-故障信号小波相干谱。正常- 故障小波相干谱中相干性较大的频带为轴承系统固有的振动信号分量频带，故障-故障小波相干谱中相干性较大的频带为轴承系统固有的振动信号分量频带和故障共振频带，由两个交叉谱对比分析可获得最优轴承故障共振频带分量。

正常-故障信号小波相干图谱

故障-故障信号小波相干图谱

clc;clear all;close all
% Signal specs
fs = 1000;                    % Sampling frequency
Flimits = [3 100];            % Frequency band of interest in Hz
nv = 12;                      % Number of logarithmic divisions per octave of frquency 
                              % (a measure of granularity of the frequency axis, usually 12/16
                              % is good enough)
alpha = 0.05;                  % significance level


%% Approach 1: This approach is similar to the Fourier-based coherence 
% Raw signals
tv = linspace(-1.5,1.5,3000); % time axis in seconds
x = rand(3000,25); % time vs. trial
y = rand(3000,25); % time vs. trial
[MsqC,SigVal,F,Coi,Phase] = WCoherence_trialAveraged(x,y,Flimits,fs,nv,alpha);


% Remove edge effects
for k = 1:size(MsqC,1)
    MsqC(k,F(k)<Coi) = NaN;
    Phase(k,F(k)<Coi) = NaN;
end


% Plot results
figure;
pcolor(tv,log2(F),MsqC); shading flat;
fvals  = [5 10 20 30 40 50 60 70 80 90];
yticks(log2(fvals));
yticklabels(fvals);
ylabel('Frequency (Hz)');
xlabel('Time (s)');
set(gca,'FontSize',20);
colormap winter; colorbar;
hold on;
ArrowDensity = [100,2]; % Density of arrows for phase indication in plot
plotPhase(gca,Phase,tv,log2(F),ArrowDensity(1),ArrowDensity(2)); 
% Create a binary matrix where values exceed the threshold
binaryMatrix = MsqC > SigVal;
contour(binaryMatrix, [1, 1], 'k'); % Contour around significant values


%% Approach 2: This approach is used clasically in wavelet coherence
% Raw signals
tv = linspace(-5,5,10000); % time axis in seconds
x = rand(1,10000); % 10 seconds of data
y = rand(1,10000); % 10 seconds of data 
[MsqC,sigvalMC,sigvalTh,F,Coi,Phase] = WCoherence_Classic(x,y,Flimits,fs,nv,10,alpha,5,6/(2*pi));


% Remove edge effects
for k = 1:size(MsqC,1)
    MsqC(k,F(k)<Coi) = NaN;
    Phase(k,F(k)<Coi) = NaN;
end


% Plot results
figure;
pcolor(tv,log2(F),MsqC); shading flat; 
fvals  = [5 10 20 30 40 50 60 70 80 90];
yticks(log2(fvals));
yticklabels(fvals);
ylabel('Frequency (Hz)');
xlabel('Time (s)');
set(gca,'FontSize',20);
colormap winter;
hold on;
ArrowDensity = [100,2]; % Density of arrows for phase indication in plot
plotPhase(gca,Phase,tv,log2(F),ArrowDensity(1),ArrowDensity(2));
% Create a binary matrix where values exceed the threshold
binaryMatrix = MsqC > sigvalMC;
contour(binaryMatrix, [1, 1], 'k'); % Contour around significant values (Monte-Carlo approach)
binaryMatrix = MsqC > sigvalTh;
%代码https://mbd.pub/o/bread/mbd-ZZmUlJdq
contour(binaryMatrix, [1, 1], 'b'); % Contour around significant values (theoretical approach)