双谱分析：振动信号的高阶统计量计算与分析

双谱分析是高阶统计量分析的重要工具，能够揭示信号中的非线性特征和相位耦合信息。

双谱分析理论基础

双谱是三阶累积量的二维傅里叶变换，能够分析信号中的二次相位耦合。与功率谱相比，双谱保留了信号的相位信息，对高斯噪声不敏感，适合分析非线性系统。

MATLAB实现代码

% 双谱分析：振动信号的高阶统计量计算与分析
clear; clc; close all;

%% 1. 生成模拟振动信号（包含非线性特征）
fs = 1000;                  % 采样频率 (Hz)
t = 0:1/fs:5;              % 时间向量 (5秒)
N = length(t);             % 信号长度

% 基本频率成分
f1 = 50;                   % 基频 (Hz)
f2 = 100;                  % 二次谐波
f3 = 150;                  % 三次谐波

% 生成信号 - 包含二次相位耦合
x = 1.0 * sin(2*pi*f1*t + 0.5) + ...
    0.6 * sin(2*pi*f2*t + 0.8) + ...
    0.4 * sin(2*pi*f3*t + 1.2) + ...
    0.3 * sin(2*pi*(f1+f2)*t + 0.5+0.8) + ...  % 相位耦合成分
    0.2 * sin(2*pi*(f1+f3)*t + 0.5+1.2);       % 相位耦合成分

% 添加高斯噪声
noise_level = 0.2;
x = x + noise_level * randn(1, N);

% 绘制时域信号
figure;
subplot(2,1,1);
plot(t, x);
xlabel('时间 (s)');
ylabel('幅值');
title('振动信号时域波形');
grid on;

% 绘制频域信号
subplot(2,1,2);
f = linspace(0, fs/2, N/2+1);
X = fft(x);
plot(f, 2*abs(X(1:N/2+1))/N);
xlabel('频率 (Hz)');
ylabel('幅值');
title('振动信号频谱');
grid on;
xlim([0, 250]);

%% 2. 计算信号的基本统计量
fprintf('===== 基本统计量 =====\n');
fprintf('信号长度: %d\n', N);
fprintf('均值: %.4f\n', mean(x));
fprintf('方差: %.4f\n', var(x));
fprintf('标准差: %.4f\n', std(x));
fprintf('偏度: %.4f\n', skewness(x));
fprintf('峰度: %.4f\n', kurtosis(x));

%% 3. 计算高阶累积量
% 三阶累积量
max_lag = 100;  % 最大滞后
[c3, lags] = bispecdx(x, max_lag);

% 绘制三阶累积量
figure;
mesh(lags, lags, c3);
xlabel('滞后 \tau_1');
ylabel('滞后 \tau_2');
zlabel('三阶累积量');
title('三阶累积量曲面');
colorbar;

%% 4. 计算双谱
% 双谱参数设置
nfft = 256;        % FFT点数
overlap = 0.5;     % 重叠率
window = hamming(nfft); % 窗函数

% 计算双谱
[bisp, waxis] = bispecd(x, nfft, window, overlap);

% 绘制双谱幅度
figure;
subplot(2,2,1);
contourf(waxis*fs/2, waxis*fs/2, abs(bisp), 20);
xlabel('频率 f1 (Hz)');
ylabel('频率 f2 (Hz)');
title('双谱幅度等高线');
colorbar;
axis square;

subplot(2,2,2);
mesh(waxis*fs/2, waxis*fs/2, abs(bisp));
xlabel('频率 f1 (Hz)');
ylabel('频率 f2 (Hz)');
zlabel('幅度');
title('双谱幅度曲面');
view(45, 30);

% 绘制双谱相位
subplot(2,2,3);
contourf(waxis*fs/2, waxis*fs/2, angle(bisp), 20);
xlabel('频率 f1 (Hz)');
ylabel('频率 f2 (Hz)');
title('双谱相位等高线');
colorbar;
axis square;

subplot(2,2,4);
mesh(waxis*fs/2, waxis*fs/2, angle(bisp));
xlabel('频率 f1 (Hz)');
ylabel('频率 f2 (Hz)');
zlabel('相位 (rad)');
title('双谱相位曲面');
view(45, 30);

%% 5. 双谱切片分析 - 检测相位耦合
% 选择特定频率切片
f1_idx = find(waxis*fs/2 >= f1, 1);  % 基频位置
f2_idx = find(waxis*fs/2 >= f2, 1);  % 二次谐波位置

% 绘制f1=f2切片
figure;
subplot(2,2,1);
slice = diag(bisp);
plot(waxis*fs/2, abs(slice));
xlabel('频率 f (Hz)');
ylabel('双谱幅度');
title('f1 = f2 切片');
grid on;

% 绘制f1=50Hz切片
subplot(2,2,2);
slice_f1 = bisp(f1_idx, :);
plot(waxis*fs/2, abs(slice_f1));
xlabel('频率 f2 (Hz)');
ylabel('双谱幅度');
title(sprintf('f1 = %.1f Hz 切片', f1));
grid on;

% 绘制f2=100Hz切片
subplot(2,2,3);
slice_f2 = bisp(:, f2_idx);
plot(waxis*fs/2, abs(slice_f2));
xlabel('频率 f1 (Hz)');
ylabel('双谱幅度');
title(sprintf('f2 = %.1f Hz 切片', f2));
grid on;

% 检测相位耦合
subplot(2,2,4);
% 寻找双谱峰值位置
[max_val, max_idx] = max(abs(bisp(:)));
[row, col] = ind2sub(size(bisp), max_idx);
f1_peak = waxis(row) * fs/2;
f2_peak = waxis(col) * fs/2;

fprintf('\n===== 双谱分析结果 =====\n');
fprintf('双谱最大峰值位置: (%.1f Hz, %.1f Hz)\n', f1_peak, f2_peak);
fprintf('峰值幅度: %.4f\n', max_val);

% 检查是否存在相位耦合
if abs(f1_peak + f2_peak - (f1+f2)) < 5 && abs(f1_peak - f1) < 5
    fprintf('检测到相位耦合: %.1f Hz + %.1f Hz = %.1f Hz\n', f1, f2, f1+f2);
elseif abs(f1_peak + f2_peak - (f1+f3)) < 5 && abs(f1_peak - f1) < 5
    fprintf('检测到相位耦合: %.1f Hz + %.1f Hz = %.1f Hz\n', f1, f3, f1+f3);
else
    fprintf('未检测到明显的相位耦合\n');
end

%% 6. 双谱积分特征提取
% 双谱幅度的积分特征
bisp_mag = abs(bisp);
total_power = sum(bisp_mag(:));
mean_mag = mean(bisp_mag(:));
std_mag = std(bisp_mag(:));

% 对角线特征
diag_mag = diag(bisp_mag);
diag_power = sum(diag_mag);
diag_mean = mean(diag_mag);

fprintf('\n===== 双谱积分特征 =====\n');
fprintf('双谱总功率: %.4f\n', total_power);
fprintf('双谱幅度均值: %.4f\n', mean_mag);
fprintf('双谱幅度标准差: %.4f\n', std_mag);
fprintf('对角线总功率: %.4f\n', diag_power);
fprintf('对角线幅度均值: %.4f\n', diag_mean);

%% 7. 双相干性分析 (Bicoherence)
% 双相干性用于量化相位耦合强度
[bic, waxis_bic] = bicoher(x, nfft, window, overlap);

figure;
subplot(1,2,1);
contourf(waxis_bic*fs/2, waxis_bic*fs/2, abs(bic), 20);
xlabel('频率 f1 (Hz)');
ylabel('频率 f2 (Hz)');
title('双相干性');
colorbar;
axis square;

subplot(1,2,2);
mesh(waxis_bic*fs/2, waxis_bic*fs/2, abs(bic));
xlabel('频率 f1 (Hz)');
ylabel('频率 f2 (Hz)');
zlabel('双相干性');
title('双相干性曲面');
view(45, 30);

%% 8. 辅助函数定义

function [c3, lags] = bispecdx(x, max_lag)
    % 计算三阶累积量
    % x: 输入信号
    % max_lag: 最大滞后
    % c3: 三阶累积量矩阵
    % lags: 滞后向量
    
    N = length(x);
    lags = -max_lag:max_lag;
    n_lags = length(lags);
    c3 = zeros(n_lags, n_lags);
    
    % 中心化信号
    x_centered = x - mean(x);
    
    % 计算三阶累积量
    for i = 1:n_lags
        tau1 = lags(i);
        for j = 1:n_lags
            tau2 = lags(j);
            
            % 确保索引不越界
            start_idx = max(1, 1 - min(0, min(tau1, tau2)));
            end_idx = min(N, N - max(0, max(tau1, tau2)));
            
            % 计算累积量
            c3(i, j) = mean(x_centered(start_idx:end_idx) .* ...
                            x_centered(start_idx+tau1:end_idx+tau1) .* ...
                            x_centered(start_idx+tau2:end_idx+tau2));
        end
    end
end

function [bisp, waxis] = bispecd(x, nfft, window, overlap)
    % 计算直接法双谱估计
    % x: 输入信号
    % nfft: FFT点数
    % window: 窗函数
    % overlap: 重叠率 (0-1)
    % bisp: 双谱矩阵
    % waxis: 频率轴
    
    N = length(x);
    window = window(:); % 确保窗函数是列向量
    L = length(window);
    
    % 计算分段数和重叠样本数
    overlap_samples = round(L * overlap);
    step = L - overlap_samples;
    n_segments = floor((N - L) / step) + 1;
    
    % 初始化双谱矩阵
    bisp = zeros(nfft, nfft);
    
    % 分段处理
    for i = 1:n_segments
        start_idx = (i-1)*step + 1;
        end_idx = start_idx + L - 1;
        
        % 提取段并加窗
        segment = x(start_idx:end_idx) .* window;
        
        % 计算FFT
        X = fft(segment, nfft);
        
        % 计算双谱估计
        for k1 = 1:nfft
            for k2 = 1:nfft
                k3 = mod(k1 + k2 - 1, nfft) + 1; % 避免索引为0
                bisp(k1, k2) = bisp(k1, k2) + X(k1) * X(k2) * conj(X(k3));
            end
        end
    end
    
    % 平均
    bisp = bisp / n_segments;
    
    % 生成频率轴
    waxis = (0:nfft/2) / (nfft/2) * 0.5;
    
    % 只保留主区域 (f1 >= 0, f2 >= 0, f1+f2 <= fs/2)
    bisp = bisp(1:nfft/2+1, 1:nfft/2+1);
end

function [bic, waxis] = bicoher(x, nfft, window, overlap)
    % 计算双相干性
    % x: 输入信号
    % nfft: FFT点数
    % window: 窗函数
    % overlap: 重叠率 (0-1)
    % bic: 双相干性矩阵
    % waxis: 频率轴
    
    N = length(x);
    window = window(:); % 确保窗函数是列向量
    L = length(window);
    
    % 计算分段数和重叠样本数
    overlap_samples = round(L * overlap);
    step = L - overlap_samples;
    n_segments = floor((N - L) / step) + 1;
    
    % 初始化矩阵
    bisp = zeros(nfft, nfft);
    P1 = zeros(nfft, 1);
    P2 = zeros(nfft, 1);
    P3 = zeros(nfft, 1);
    
    % 分段处理
    for i = 1:n_segments
        start_idx = (i-1)*step + 1;
        end_idx = start_idx + L - 1;
        
        % 提取段并加窗
        segment = x(start_idx:end_idx) .* window;
        
        % 计算FFT
        X = fft(segment, nfft);
        P = abs(X).^2;
        
        % 累积功率谱
        P1 = P1 + P;
        
        % 计算双谱估计
        for k1 = 1:nfft
            for k2 = 1:nfft
                k3 = mod(k1 + k2 - 1, nfft) + 1; % 避免索引为0
                bisp(k1, k2) = bisp(k1, k2) + X(k1) * X(k2) * conj(X(k3));
                
                % 累积功率乘积
                if i == 1
                    P2(k1) = P2(k1) + P(k1) * P(k2);
                    P3(k1) = P3(k1) + P(k1) * P(k3);
                end
            end
        end
    end
    
    % 平均
    bisp = bisp / n_segments;
    P1 = P1 / n_segments;
    P2 = P2 / n_segments;
    P3 = P3 / n_segments;
    
    % 计算双相干性
    bic = zeros(nfft, nfft);
    for k1 = 1:nfft
        for k2 = 1:nfft
            k3 = mod(k1 + k2 - 1, nfft) + 1;
            denominator = sqrt(P1(k1) * P1(k2) * P1(k3));
            if denominator > 0
                bic(k1, k2) = abs(bisp(k1, k2)) / denominator;
            else
                bic(k1, k2) = 0;
            end
        end
    end
    
    % 生成频率轴
    waxis = (0:nfft/2) / (nfft/2) * 0.5;
    
    % 只保留主区域 (f1 >= 0, f2 >= 0, f1+f2 <= fs/2)
    bic = bic(1:nfft/2+1, 1:nfft/2+1);
end

function [bisp, waxis] = bispecd_indirect(x, nfft, max_lag)
    % 间接法双谱估计（通过三阶累积量）
    % x: 输入信号
    % nfft: FFT点数
    % max_lag: 最大滞后
    % bisp: 双谱矩阵
    % waxis: 频率轴
    
    % 计算三阶累积量
    [c3, lags] = bispecdx(x, max_lag);
    
    % 对三阶累积量进行二维FFT
    bisp = fft2(c3, nfft, nfft);
    
    % 生成频率轴
    waxis = (0:nfft/2) / (nfft/2) * 0.5;
    
    % 只保留主区域 (f1 >= 0, f2 >= 0, f1+f2 <= fs/2)
    bisp = bisp(1:nfft/2+1, 1:nfft/2+1);
end

代码说明

这个MATLAB实现提供了完整的双谱分析流程，包括：

1. 信号生成：创建包含相位耦合成分的模拟振动信号
2. 基本统计量计算：均值、方差、偏度、峰度等
3. 三阶累积量计算：双谱分析的基础
4. 双谱估计：使用直接法和间接法计算双谱
5. 双谱可视化：幅度和相位的等高线图和曲面图
6. 切片分析：检测特定频率的相位耦合
7. 双相干性分析：量化相位耦合强度
8. 特征提取：从双谱中提取有用的特征

关键功能

1. 双谱估计：

   • bispecd函数使用直接法计算双谱
   • bispecd_indirect函数使用间接法（通过三阶累积量）计算双谱

2. 双相干性计算：

   • bicoher函数计算双相干性，用于量化相位耦合强度

3. 三阶累积量计算：

   • bispecdx函数计算信号的三阶累积量

参考代码 双谱分析，对振动信号进行高阶统计量的计算分析 www.youwenfan.com/contentcnj/63908.html

应用

双谱分析在以下领域有广泛应用：

1. 机械故障诊断：检测旋转机械中的非线性特征和故障
2. 生物医学信号处理：分析EEG、ECG等生物信号中的相位耦合
3. 地球物理信号处理：分析地震信号中的非线性特征
4. 通信系统：分析非线性失真和相位耦合

使用

1. 代码首先生成一个包含相位耦合的模拟振动信号
2. 然后计算并显示信号的基本统计量
3. 接着计算三阶累积量和双谱，并进行可视化
4. 通过切片分析检测相位耦合
5. 最后计算双相干性并提取双谱特征

注意

1. 双谱计算量较大，对于长信号可能需要较长时间
2. 参数选择（如FFT点数、窗函数、重叠率）会影响分析结果
3. 双谱分析对非高斯和非线性信号最有效
4. 实际应用中可能需要预处理信号（如去趋势、滤波等）

这个实现提供了一个完整的双谱分析框架，您可以根据实际振动信号分析需求进行修改和扩展。