%%%%%%%%% 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)]); % 注释掉这部分