%%%%%%%%% the signal with fast oscillating IF %%%%%%%%%%%%%
% Chen S, Yang Y, Peng Z, et al, Detection of Rub-Impact Fault for Rotor-Stator Systems: A Novel Method Based on Adaptive Chirp Mode Decomposition, Journal of Sound and Vibration, 2018.
% instantaneous frequency initialization by finding peaks of Fourier spectrum
clc
clear
close all
% 提示用户输入信号文件的路径
filename = input('请输入信号文件的完整路径:', 's');
% 检查文件是否存在
if exist(filename, 'file')
% 文件存在,读取信号
Sig = load(filename);
if length(size(Sig)) > 1
Sig = Sig(:, 1); % 假设信号在第一列
end
% 根据信号的实际长度重新生成时间向量
t = linspace(0, 1, length(Sig));
else
error('文件不存在,请检查路径是否正确。');
end
% 绘制原始信号
figure
set(gcf,'Position',[20 100 640 500]);
set(gcf,'Color','w');
plot(t,Sig,'linewidth',2);
xlabel('Time / Sec','FontSize',24,'FontName','Times New Roman');
ylabel('Amplitude / AU','FontSize',24,'FontName','Times New Roman');
set(gca,'FontSize',24)
set(gca,'linewidth',2);
%% Fourier spectrum
Spec = 2abs(fft(Sig))/length(Sig);
Spec = Spec(1:round(end/2));
Freqbin = linspace(0, 1000/2, length(Spec));
figure
set(gcf,'Position',[20 100 640 500]);
set(gcf,'Color','w');
plot(Freqbin, Spec, 'linewidth', 2);
xlabel('Frequency / Hz','FontSize',24,'FontName','Times New Roman');
ylabel('Amplitude / AU','FontSize',24,'FontName','Times New Roman');
set(gca,'FontSize',24)
set(gca,'linewidth',2);
axis([100 500 0 1.1max(Spec)])
%% STFT
% 设置 STFT 参数
window = hamming(min(512, length(Sig))); % 使用较小的窗口大小
noverlap = round(min(256, length(Sig) / 2)); % 确保重叠不是信号的一半以上
nfft = min(512, length(Sig)); % 使用较小的 FFT 长度
% 调用 STFT 函数
[TFSpec, f, t_stft] = stft(Sig, 'Window', window, 'OverlapLength', noverlap, 'FFTLength', nfft);
% 绘制 STFT 结果
figure
imagesc(t_stft, f, abs(TFSpec));
axis xy;
axis([0 max(t_stft) 0 max(f)]);
set(gcf,'Position',[20 100 320 250]);
xlabel('Time / Sec','FontSize',12,'FontName','Times New Roman');
ylabel('Frequency / Hz','FontSize',12,'FontName','Times New Roman');
set(gca,'YDir','normal')
set(gca,'FontSize',12);
set(gcf,'Color','w');
%% parameter setting
alpha0 = 1e-3;
beta = 1e-4;
tol = 1e-8;
%% Component 1 extraction
[~, findex1] = max(Spec);
peakfre1 = Freqbin(findex1);
iniIF1 = peakfre1 * ones(1, length(Sig));
[IFest1, IAest1, sest1] = ACMD(Sig, 1000, iniIF1, alpha0, beta, tol);
%% estimated IF
% 确保 IFest1 与时间向量 t 的长度一致
IFest1 = IFest1(1:length(t));
figure
plot(t, IFest1, 'r--', 'linewidth', 3);
set(gcf,'Position',[20 100 640 500]);
xlabel('Time / Sec','FontSize',24,'FontName','Times New Roman');
ylabel('Frequency / Hz','FontSize',24,'FontName','Times New Roman');
set(gca,'YDir','normal')
set(gca,'FontSize',24);
set(gca,'linewidth',2);
set(gca,'FontName','Times New Roman')
set(gcf,'Color','w');
axis([0 1 100 500]);
%% Reconstructed modes
figure
set(gcf,'Position',[20 100 640 500]);
subplot(2,1,1)
set(gcf,'Color','w');
plot(t,sest1,'b--','linewidth',2);
ylabel('C1','FontSize',24,'FontName','Times New Roman');
set(gca,'YDir','normal')
set(gca,'FontSize',24);
set(gca,'linewidth',2);
axis([0 1 -1 1])
% 注释掉 Component 2 的提取和绘图
% subplot(2,1,2)
% set(gcf,'Color','w');
% plot(t,sest2,'b--','linewidth',2);
% ylabel('C2','FontSize',24,'FontName','Times New Roman');
% xlabel('Time / Sec','FontSize',24,'FontName','Times New Roman')
% set(gca,'YDir','normal')
% set(gca,'FontSize',24);
% set(gca,'linewidth',2);
% axis([0 1 -1 1])
%% adaptive time-frequency spectrum
band = [0 500];
IFmulti = [IFest1]; % 只使用 IFest1
IAmulti = [IAest1]; % 只使用 IAest1
[ASpec, fbin] = TFspec(IFmulti, IAmulti, band);
figure
imagesc(t, fbin, abs(ASpec));
axis([0 1 100 500]);
set(gcf,'Position',[20 100 640 500]);
xlabel('Time / Sec','FontSize',24,'FontName','Times New Roman');
ylabel('Frequency / Hz','FontSize',24,'FontName','Times New Roman');
set(gca,'YDir','normal')
set(gca,'FontSize',24);
set(gcf,'Color','w');
colorbar('YTickLabel',{' '})
% 计算并显示SNR
SNR1 = SNR(Sig, sest1);
% SNR2 = SNR(Sig, sest2); % 注释掉这部分
disp(['SNR for Component 1: ', num2str(SNR1)]);
% disp(['SNR for Component 2: ', num2str(SNR2)]); % 注释掉这部分
在处理地震动数据并绘制图像时,如果您发现绘制的图像全部接近于零,这可能是由于多个因素导致的。以下是一些可能的问题和建议,您可以根据这些信息进行调节和排查。
首先,确保您加载的数据在合理的范围内。您可以通过打印信号的最小值和最大值来检查:
disp(['Min value: ', num2str(min(Sig))]);
disp(['Max value: ', num2str(max(Sig))]);
如果信号的幅度非常小(接近于零),则可能需要考虑信号的放大。
如果信号中存在噪声,可能会影响后续的处理。考虑对信号进行滤波处理,例如使用低通滤波器来去除高频噪声。
在计算 FFT 时,确保您正确处理了 FFT 的结果。您在计算频谱时使用了以下代码:
Spec = 2 * abs(fft(Sig)) / length(Sig);
确保您取的是前半部分的频谱,并且在绘图时设置的 y 轴范围是否合理。
在绘制频谱时,您设定的 y 轴范围可能不合适:
axis([100 500 0 1.1 * max(Spec)]);
如果 max(Spec)
的值非常小,那么绘图的 y 轴范围可能会导致图像几乎不可见。可以尝试将 y 轴范围设置为更大的值,例如:
axis([100 500 0 0.1]); % 根据实际情况调整
在进行短时傅里叶变换 (STFT) 时,确保窗口大小和重叠设置合理:
window = hamming(min(512, length(Sig)));
noverlap = round(min(256, length(Sig) / 2));
nfft = min(512, length(Sig));
如果信号长度较短,可能导致 STFT 结果不明显。尝试增大窗口大小或调整重叠参数。
在调用 ACMD
函数时,您使用了以下参数:
alpha0 = 1e-3;
beta = 1e-4;
tol = 1e-8;
这些参数可能会影响到分解的结果。可以尝试调整这些参数,尤其是 alpha0
和 beta
,以观察对结果的影响。通常,较大的 alpha0
和 beta
值可能会导致更明显的模式提取。
在计算 SNR 时,确保 sest1
和 Sig
的维度匹配。如果 sest1
过小或接近零,SNR 计算的结果可能会不准确。
在调试时,可以逐步查看每个步骤的输出,例如:
disp
或 fprintf
打印中间结果,帮助定位问题。这里是一个简化的示例,帮助您进行调试:
% 检查信号
disp(['Min value: ', num2str(min(Sig))]);
disp(['Max value: ', num2str(max(Sig))]);
% 绘制频谱
Spec = 2 * abs(fft(Sig)) / length(Sig);
Spec = Spec(1:round(end/2));
Freqbin = linspace(0, 1000/2, length(Spec));
figure;
plot(Freqbin, Spec, 'linewidth', 2);
xlabel('Frequency / Hz');
ylabel('Amplitude / AU');
axis([100 500 0 0.1]); % 根据实际情况调整
% 计算和绘制 STFT
[TFSpec, f, t_stft] = stft(Sig, 'Window', window, 'OverlapLength', noverlap, 'FFTLength', nfft);
figure;
imagesc(t_stft, f, abs(TFSpec));
axis xy;
xlabel('Time / Sec');
ylabel('Frequency / Hz');
通过以上步骤和建议,您可以逐步排查导致图像接近于零的原因,并调整相关参数以获得更好的结果。如果问题仍然存在,可以提供更多的上下文或数据,以便进行进一步的分析和建议。