将时间序列转成图像——小波变换方法 Matlab实现

Wf(a,b)=<f,ψa,b>=1|a|∫−∞+∞f(t)ψ(t−ba)dt" role="presentation">Wf(a,b)=<f,ψa,b>=1|a|√∫+∞−∞f(t)ψ(t−ba)dt$W_f(a,b)=<f,{\psi }_{a,b}>=\frac{1}{\sqrt{|a|}}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}f(t)\psi \left(\frac{t-b}{a}\right)dt$​

clear, close all
 
%% initialize parameters
samplerate=500; % in Hz
 
fstep=1;   % frequency step for wavelet
 
%% generate simulated signals with step changes in frequency
data = csvread('3_1_link6_28_5_30min.csv');     %  input the signal from the Excle
 
data = data';   % change the signal from column to row
 
N = length(data);   % calculate the length of the data
 
taxis = [1:N]/samplerate;   % time axis for whole data length
 
figure, 
plot(taxis,data),xlim([taxis(1) taxis(end)])
xlabel('Time (s)')
 
%% Time-frequency analysis (CWT, morlet wavelet)
spec = tfa_morlet(data, samplerate, 1, 250, fstep);
faxis=[1:fstep:250];
Mag=abs(spec);     % get spectrum magnitude
 
im = figure('color',[1 1 1]);
imagesc(taxis,faxis,Mag)   % plot spectrogram as an image
colorbar
axis([taxis(1) taxis(end) faxis(1) faxis(end)])
xlabel('Time (s)')
ylabel('Frequency (Hz)')
title('Time-frequency analysis (CWT)')
 
saveas(im,'CWT_1.bmp')
 
function TFmap = tfa_morlet(td, fs, fmin, fmax, fstep)
TFmap = [];
for fc=fmin:fstep:fmax
    MW = MorletWavelet(fc/fs);  % calculate the Morlet Wavelet by giving the central freqency
    cr = conv(td, MW, 'same');  % convolution
    
    TFmap = [TFmap; abs(cr)];
end
 
function MW = MorletWavelet(fc)
 
F_RATIO = 7;    % frequency ratio (number of cycles): fc/sigma_f, should be greater than 5
Zalpha2 = 3.3;  % value of Z_alpha/2, when alpha=0.001
 
sigma_f = fc/F_RATIO;
sigma_t = 1/(2*pi*sigma_f);
A = 1/sqrt(sigma_t*sqrt(pi));
max_t = ceil(Zalpha2 * sigma_t);
 
t = -max_t:max_t;
 
%MW = A * exp((-t.^2)/(2*sigma_t^2)) .* exp(2i*pi*fc*t);
v1 = 1/(-2*sigma_t^2);
v2 = 2i*pi*fc;
MW = A * exp(t.*(t.*v1+v2));