小波变换：多分辨率的时频分析

% 亲手画一个 10Hz 的 Morlet 小波
fs = 500; f0 = 10; n_cyc = 6;
sigma_t = n_cyc / (2*pi*f0);          % 高斯包络的标准差由周期数决定
t = -1 : 1/fs : 1;
w = exp(2i*pi*f0*t) .* exp(-t.^2/(2*sigma_t^2));

figure;
subplot(2,1,1); plot(t, real(w), 'b', t, abs(w), 'r--');
legend('实部（振荡）','包络'); title('Morlet 小波 = 高斯包络 × 正弦');
xlim([-0.4 0.4]);
subplot(2,1,2); plot(t, abs(w), 'r'); title('包络：决定有效窗长');
xlim([-0.4 0.4]);

% Morlet：频率 ≈ 中心频率 / 尺度（以母小波中心频率 f_c 归一）
% Matlab cwt 默认 Morse 小波时可用scal2freq换算：
%   fb = 3,  Morse 近似 Morlet（n_cyc 大致对应）
f_c = 1;                          % 母小波中心频率（示例）
scales = [0.05 0.1 0.25 0.5 1 2 4 8];
freqs  = f_c ./ scales;           % 频率 = 中心频率 / 尺度
disp(table(scales', freqs', 'VariableNames', {'scale','freq_Hz'}));
% 直觉：尺度是"拉长倍数"。拉长 2 倍 → 周期变 2 倍 → 频率减半

% x : 一维信号；fs : 采样率
fs = 500;
[cfs, frq] = cwt(x, 'morse3', fs);     % Morse(3) 近似经典 Morlet
figure;
surface(seconds(0:1/fs:(numel(x)-1)/fs), frq, abs(cfs));
axis xy; ylim([2 45]);
xlabel('时间 (s)'); ylabel('频率 (Hz)');
title('连续小波变换幅度谱');
colormap(parula); colorbar;

% 也可以直接用便捷绘图接口：
% cwt(x, fs)  自动配色与频率轴

function [P, phase, f, t] = morletTFR(x, fs, fList, n_cyc)
% x: 1 x N 信号；fList: 频率列表（如 2:1:45）；n_cyc: 周期数
x = x(:).';  N = numel(x);
t = (0:N-1) / fs;
nF = numel(fList);  P = zeros(nF, N);  phase = P;

X = fft(x);                                  % 信号只变换一次
for k = 1:nF
    f0  = fList(k);
    s_t = n_cyc / (2*pi*f0);                 % 该频率下的包络宽度
    tw  = -3*s_t : 1/fs : 3*s_t;             % 覆盖 ±3σ
    w   = exp(2i*pi*f0*tw) .* exp(-tw.^2/(2*s_t^2));
    w   = w / sqrt(s_t);                     % 能量归一化

    W = fft(w, N);                           % 小波补零到信号长度
    convr = ifft( X .* W );                   % 频域相乘 = 时域卷积
    convr = [convr(end-floor(numel(w)/2)+1:end), convr(1:end-floor(numel(w)/2))];
    P(k,:)     = abs(convr).^2;              % 瞬时功率
    phase(k,:) = angle(convr);               % 瞬时相位（8.7 ITPC 会用到）
end
end

% 用法示例：
% [P, ph, f, t] = morletTFR(eeg_chan, 500, 2:1:45, 6);
% imagesc(t, f, 10*log10(P)); axis xy;