目录
数据
首先我们需要准备一个mat数据文件。文件格式如下
这个项目提供了24中转图像的类型：
时频类：
梅尔频谱图Mel spectrogram
短时傅里叶变换short-time Fourier transform
变换S-transform
魏格纳分布Wigner-Ville Distribution
离散魏格纳分布Discrete Wigner-Ville Distribution
希尔伯特变换Hilbert-Huang Transform
连续小波变换Continuous wavelet transform
实小波变换Real wavelet transform
同步压缩变换Synchrosqueezing transform
小波同步压缩变换wavelet synchrosqueezed transform
小波二阶同步压缩变换wavelet second order synchrosqueezed transform
垂直二阶同步压缩变换vertical second-order synchrosqueezing
多尺度同步压缩变换Multisynchrosqueezing Transform
小波多尺度同步压缩变换Wavelet Multisynchrosqueezed Transform
局部最大同步压缩变换Local maximum synchrosqueezing transform
时间重分配多同步压缩变换Time-reassigned Multisynchrosqueezing Transform
同步提取变换Synchroextracted transform
小波同步提取变换Wavelet Synchroextracted Transform
暂态提取变换transient-extracting transform
二阶暂态提取变换Second-order transient-extracting transform
转换类：
格拉姆角和场Gramian angular summation field
递归图recurrence plots
相对位置矩阵Relative Position Matrix
下面是格拉姆角差场的样例图
递归图
变换S
我们要创建以下文件。
main.m
%% 1D signal to 2D image 实验
clc; clear; close all
%% 1.加载数据 → 一维数据
load test_data
%% 2. 选择方法
n_point = 1024; % 以1024 个点 划分样本
flag =0; % 坐标图窗控制参数，flag=0 时显示坐标信息，=1 不显示坐标信息
method = 'ST'; % 字符串格式
% method可选：
% 时频类：
% method = 'mel'， 运行→ 梅尔频谱图Mel spectrogram
% method = 'stft'， 运行→ 短时傅里叶变换short-time Fourier transform
% method = 'ST'， 运行→ s变换S-transform
% method = 'WVD'， 运行→ 魏格纳分布Wigner-Ville Distribution
% method = 'DWVD'， 运行→ 离散魏格纳分布Discrete Wigner-Ville Distribution
% method = 'HHT'， 运行→ 希尔伯特变换Hilbert-Huang Transform
% method = 'cwt'， 运行→ 连续小波变换Continuous wavelet transform
% method = 'RWT'， 运行→ 实小波变换Real wavelet transform
% method = 'SST'， 运行→ 同步压缩变换Synchrosqueezing transform
% method = 'wsst'， 运行→ 小波同步压缩变换wavelet synchrosqueezed transform
% method = 'wsst2'， 运行→ 小波二阶同步压缩变换wavelet second order synchrosqueezed transform
% method = 'VSST2'， 运行→ 垂直二阶同步压缩变换vertical second-order synchrosqueezing
% method = 'MSST'， 运行→ 多尺度同步压缩变换Multisynchrosqueezing Transform
% method = 'wmsst'， 运行→ 小波多尺度同步压缩变换Wavelet Multisynchrosqueezed Transform
% method = 'LMSST'， 运行→ 局部最大同步压缩变换Local maximum synchrosqueezing transform
% method = 'TMSST'， 运行→ 时间重分配多同步压缩变换Time-reassigned Multisynchrosqueezing Transform
% method = 'SET'， 运行→ 同步提取变换Synchroextracted transform
% method = 'wset'， 运行→ 小波同步提取变换Wavelet Synchroextracted Transform
% method = 'TET'， 运行→ 暂态提取变换transient-extracting transform
% method = 'STET'， 运行→ 二阶暂态提取变换Second-order transient-extracting transform
% 转换类：
% method = 'GASF'， 运行→ 格拉姆角和场Gramian angular summation field
% method = 'GADF'， 运行→ 格拉姆角差场Gramian angular difference field
% method = 'RP'， 运行→ 递归图recurrence plots
% method = 'RPM'， 运行→ 相对位置矩阵Relative Position Matrix
%% 3. 定义文件路径，用于保存图像
%% 这里是自动定义的文件路径，如不需要可将此部分注释后，手动添加路径，类似 file_path='D:\code\...'
%% file_path 一定要给文件路径
file_path =[method '_images']; %
% 判断文件路径是否存在，如果文件路径不存在，则创建
if ~exist(file_path, 'dir')
mkdir(file_path);
fprintf('文件路径 %s 创建成功！\n', file_path);
end
%% 4.生成图像
addpath(genpath('Algorithms')) % 将算法加入路径
get_images(x,n_point,method,file_path,flag);
rmpath(genpath('Algorithms')) % 使用完后，将算法移除路径1.加载数据
load test_data
2.选择方法
method = 'ST';
3.文件路径管理
file_path =[method '_images'];
if ~exist(file_path, 'dir')
mkdir(file_path);
end
自动创建存储转换后 2D 图像的目录，例如 ST_images/。
4.生成图像
addpath(genpath('Algorithms')) % 加载算法
get_images(x,n_point,method,file_path,flag); % 生成并保存图像
rmpath(genpath('Algorithms')) % 移除算法路径
关键函数 get_images(x,n_point,method,file_path,flag) 负责信号处理和图像生成
get_images.m
这段 MATLAB 代码是一个函数 get_images，用于对输入信号 x 进行不同的时频分析方法（如 Mel 频谱、短时傅里叶变换 STFT、小波变换 CWT、同步压缩变换 SST 等），并将生成的时频图像保存为 jpg 文件。
%% 注意：以下参数只是测试用，需结合自身专业知识设定
function get_images(x,n_point,method,file_path,flag)
% 包括步骤：调用算法，生成图形，保存
x=x(:); % 统一转换成 一列 数据
n_sample = floor(length(x)/n_point); % 样本个数， floor是为了防止原始信号长度不能整除 n_point
switch method
case 'mel'
% 参数设置
fs = 5000; % 采样频率
WL=10; % 一次分析选取的信号长度, <输入信号长度
OL=5; % 分析窗口重叠长度，即两个相邻窗口的重叠长度 < WL
FL=1024; % 计算DFT的点数，大于或等于 WL 的正整数
for i=1:n_sample
xx = x(1+n_point*(i-1):i*n_point,:); % 循环读取每个样本
[S,cF,t] =melSpectrogram(xx,fs,'WindowLength',WL,'OverlapLength',OL,'FFTLength',FL); % S为计算得到的mel频谱图
S = 10*log10(S+eps); % 转换为dB
% 绘图
imagesc(t,cF,S)
axis_control(flag) % 坐标轴控制
% 保存
saveas(gcf,fullfile(file_path, [method '_' num2str(i) '.jpg'])) % 保存
end
case 'stft'
% STFT参数
fs=5000; %采样频率
window_len=10; %设置窗口长度。
window=hamming(window_len);%设置汉宁窗，当然也可以使用其他的窗
overlap = 5; % 重叠数，overlap<window_len
nfft = 1024; % DFT 点数
% 计算 STFT
for i=1:n_sample
xx = x(1+n_point*(i-1):i*n_point,:); % 循环读取每个样本
[S,F,T,~]=spectrogram(xx,window,overlap,nfft,fs);
%绘图
imagesc(T, F, 20*log10(abs(S)));
axis_control(flag)
% 保存
saveas(gcf,fullfile(file_path, [method '_' num2str(i) '.jpg'])) % 保存
end
case 'cwt'
% 小波参数
fs=5000; % 采样频率
wavename='haar'; % 选用haar，可以替换其他的小波函数
totalscal=256; %尺度序列的长度，即scal的长度
% 其他参数 不用修改
t=1/fs:1/fs:1;
fc=centfrq(wavename); %小波的中心频率
cparam=2*fc*totalscal;
a=totalscal:-1:1;
scal=cparam./a;
% 计算 cwt
for i=1:n_sample
xx = x(1+n_point*(i-1):i*n_point,:); % 循环读取每个样本
coefs=cwt(xx,scal,wavename); %得到小波系数
f=scal2frq(scal,wavename,1/fs); %将尺度转换为频率
coef=abs(coefs);
%绘图
imagesc(t,f,coef);
axis_control(flag)
% 保存
saveas(gcf,fullfile(file_path, [method '_' num2str(i) '.jpg'])) % 保存
end
case 'ST'
% s变换
for i=1:n_sample
xx = x(1+n_point*(i-1):i*n_point,:); % 循环读取每个样本
[st_matrix,st_t,st_freq] = st(xx);% 调用 st.m函数 计算
%绘图
imagesc(st_t,st_freq,abs(st_matrix));%绘制色谱图
axis_control(flag)
% 保存
saveas(gcf,fullfile(file_path, [method '_' num2str(i) '.jpg'])) % 保存
end
case 'HHT'
% HHT变换
Fs=1024;%采样频率
for i=1:n_sample
xx = x(1+n_point*(i-1):i*n_point,:); % 循环读取每个样本
% 对信号进行EMD分解
imf = emd(xx);
% hht绘制希尔伯特谱
hht(imf,Fs);
axis_control(flag)
% 保存
saveas(gcf,fullfile(file_path, [method '_' num2str(i) '.jpg'])) % 保存
end
case 'GASF'
% 格拉姆角和场
for i=1:n_sample
xx = x(1+n_point*(i-1):i*n_point,:); % 循环读取每个样本
% 计算
GM = GAF(xx,'+',0);
%绘制
imagesc(GM);
axis_control(flag)
% 保存
saveas(gcf,fullfile(file_path, [method '_' num2str(i) '.jpg'])) % 保存
end
case 'GADF'
% 格拉姆角差场
for i=1:n_sample
xx = x(1+n_point*(i-1):i*n_point,:); % 循环读取每个样本
% 计算
GM = GAF(xx,'-',0);
%绘制
imagesc(GM);
axis_control(flag)
% 保存
saveas(gcf,fullfile(file_path, [method '_' num2str(i) '.jpg'])) % 保存
end
case 'RP'
% 递归图
% 参数设置
m=4; % 嵌入维数
tau=1; % 时延
eta = 0.1; % 阈值
for i=1:n_sample
xx = x(1+n_point*(i-1):i*n_point,:); % 循环读取每个样本
% 计算
rp_mat = RecurrencePlot(xx,m,tau,eta,0);
%绘制
imshow(rp_mat,[]);%黑白
axis_control(flag)
% 保存
saveas(gcf,fullfile(file_path, [method '_' num2str(i) '.jpg'])) % 保存
end
case 'SST'
% Synchrosqueezing transform 同步压缩变换
fs = 100; % 采样频率
for i=1:n_sample
xx = x(1+n_point*(i-1):i*n_point,:); % 循环读取每个样本
time=(1:length(xx))/fs;
fre=(fs/2)/(length(xx)/2):(fs/2)/(length(xx)/2):(fs/2);
% 计算
Ts = SST(xx); % 调用 SST.m 计算
%绘制
imagesc(time,fre,abs(Ts));
axis_control(flag)
% 保存
saveas(gcf,fullfile(file_path, [method '_' num2str(i) '.jpg'])) % 保存
end
case 'SET'
% Synchroextracted transform 同步提取变换
fs = 100; % 采样频率
time=(1:n_point)/fs;
fre=(fs/2)/(n_point/2):(fs/2)/(n_point/2):(fs/2);
for i=1:n_sample
xx = x(1+n_point*(i-1):i*n_point,:); % 循环读取每个样本
% 计算
Ts = SET(xx); % 调用 SET.m 计算
%绘制
imagesc(time,fre,abs(Ts));
axis_control(flag)
% 保存
saveas(gcf,fullfile(file_path, [method '_' num2str(i) '.jpg'])) % 保存
end
case 'RPM'
% 相对位置矩阵(Relative Position Matrix)
k=4;
for i=1:n_sample
xx = x(1+n_point*(i-1):i*n_point,:); % 循环读取每个样本
% 计算
RPM = RelativePositionMatrix(xx,k); % 调用 RelativePositionMatrix.m 计算
%绘制
imagesc(RPM);
axis_control(flag)
% 保存
saveas(gcf,fullfile(file_path, [method '_' num2str(i) '.jpg'])) % 保存
end
case 'wsst'
% WSST : the synchrosqueezed transform 基于小波同步压缩变换
for i=1:n_sample
xx = x(1+n_point*(i-1):i*n_point,:); % 循环读取每个样本
% 计算-
WSST = wsst(xx);
imagesc(abs(WSST));
axis_control(flag)
% 保存
saveas(gcf,fullfile(file_path, [method '_' num2str(i) '.jpg'])) % 保存
end
case 'wsst2'
% WSST2: the second order % 基于小波二阶同步压缩变换
for i=1:n_sample
xx = x(1+n_point*(i-1):i*n_point,:); % 循环读取每个样本
% 计算-调用 Wsst2_new.m 计算
[WSST2, fs, as, ~, ~, ~, ~, ~] = Wsst2_new(xx);
imagesc(as,fs,abs(WSST2));
axis_control(flag)
% 保存
saveas(gcf,fullfile(file_path, [method '_' num2str(i) '.jpg'])) % 保存
end
case 'VSST2'
% VSST2: vertical second-order synchrosqueezing 垂直二阶同步压缩变换
% 参数设置
for i=1:n_sample
xx = x(1+n_point*(i-1):i*n_point,:); % 循环读取每个样本
% 计算-调用 VSST2.m 计算
[VSST,~,~,~,~] = VSST2(xx);
% 绘图
imagesc(abs(VSST));
axis_control(flag)
% 保存
saveas(gcf,fullfile(file_path, [method '_' num2str(i) '.jpg'])) % 保存
end
case 'STET'
% Second-order transient-extracting transform 二阶暂态提取变换
hlength=128;
for i=1:n_sample
xx = x(1+n_point*(i-1):i*n_point,:); % 循环读取每个样本
% 计算
[Te] = STET(xx,hlength);
% 绘图
imagesc(abs(Te));
axis_control(flag)
% 保存
saveas(gcf,fullfile(file_path, [method '_' num2str(i) '.jpg'])) % 保存
end
case 'TET'
% transient-extracting transform 暂态提取变换
hlength=128;
for i=1:n_sample
xx = x(1+n_point*(i-1):i*n_point,:); % 循环读取每个样本
% 计算
[Te] = TET(xx,hlength);
% 绘图
imagesc(abs(Te));
axis_control(flag)
% 保存
saveas(gcf,fullfile(file_path, [method '_' num2str(i) '.jpg'])) % 保存
end
case 'wset'
% Wavelet Synchroextracted Transform 小波同步提取变换
for i=1:n_sample
xx = x(1+n_point*(i-1):i*n_point,:); % 循环读取每个样本
% 计算
[sst,f] = wset(xx);
% 绘图
imagesc(abs(sst));
axis_control(flag)
% 保存
saveas(gcf,fullfile(file_path, [method '_' num2str(i) '.jpg'])) % 保存
end
case 'wmsst'
% Wavelet Multisynchrosqueezed Transform 小波多尺度同步变换
num=1; %1-4
for i=1:n_sample
xx = x(1+n_point*(i-1):i*n_point,:); % 循环读取每个样本
% 计算
[sst,f] = wmsst(xx,num);
% 绘图
imagesc(abs(sst));
axis_control(flag)
% 保存
saveas(gcf,fullfile(file_path, [method '_' num2str(i) '.jpg'])) % 保存
end
case 'LMSST'
% Local maximum synchrosqueezing transform 局部最大同步压缩变换
hlength=64;
for i=1:n_sample
xx = x(1+n_point*(i-1):i*n_point,:); % 循环读取每个样本
% 计算
[Ts, IF, tfr] = LMSST(xx,hlength);
% 绘图
imagesc(abs(Ts));
axis_control(flag)
% 保存
saveas(gcf,fullfile(file_path, [method '_' num2str(i) '.jpg'])) % 保存
end
case 'MSST'
% Multisynchrosqueezi?ng Transform 多同步压缩变换
for i=1:n_sample
xx = x(1+n_point*(i-1):i*n_point,:); % 循环读取每个样本
% 计算
[Ts,~,~] = MSST_new(xx);
% 绘图
imagesc(abs(Ts));
axis_control(flag)
% 保存
saveas(gcf,fullfile(file_path, [method '_' num2str(i) '.jpg'])) % 保存
end
case 'TMSST'
% Time-reassigned Multisynchrosqueezing Transform 时间重分配多同步压缩变换
for i=1:n_sample
xx = x(1+n_point*(i-1):i*n_point,:); % 循环读取每个样本
% 计算
[Ts,~] = TMSST(xx);
% 绘图
imagesc(abs(Ts));
axis_control(flag)
% 保存
saveas(gcf,fullfile(file_path, [method '_' num2str(i) '.jpg'])) % 保存
end
case 'RWT'
% Real Wavelet Transform 实小波变换
nvoice=8;
wavelet_name = 'Morlet'; % string 'Gauss', 'DerGauss','Sombrero', 'Morlet'
for i=1:n_sample
xx = x(1+n_point*(i-1):i*n_point,:); % 循环读取每个样本
% 计算-调用RWT.m计算
Ts = RWT(xx,nvoice,wavelet_name);
% 绘图
imagesc(abs(Ts));
axis_control(flag)
% 保存
saveas(gcf,fullfile(file_path, [method '_' num2str(i) '.jpg'])) % 保存
end
case 'WVD'
% Wigner-Ville Distribution 魏格纳分布
for i=1:n_sample
xx = x(1+n_point*(i-1):i*n_point,:); % 循环读取每个样本
% 计算
Ts = wvd(xx);
% 绘图
imagesc(abs(Ts));
axis_control(flag)
% 保存
saveas(gcf,fullfile(file_path, [method '_' num2str(i) '.jpg'])) % 保存
end
case 'DWVD'
% Discrete Wigner-Ville Distribution 离散魏格纳分布
fs = 1000; % 采样率
for i=1:n_sample
xx = x(1+n_point*(i-1):i*n_point,:); % 循环读取每个样本
% 计算-调用WignerDist.m计算
[~,DVW,~,~] = DWVD(xx,fs);
% 绘图
imagesc(abs(DVW));
axis_control(flag)
% 保存
saveas(gcf,fullfile(file_path, [method '_' num2str(i) '.jpg'])) % 保存
end
end
close % 关闭图窗
end1. 函数输入参数
function get_images(x, n_point, method, file_path, flag)
x：输入信号（一维数据）。n_point：每个样本的点数，即对x进行切片的大小。method：选择时频分析方法（如mel、stft、cwt等）。file_path：保存图像的路径。flag：用于控制坐标轴显示方式（由axis_control(flag)控制）。
2. 数据预处理
x = x(:); % 确保输入信号为列向量
n_sample = floor(length(x) / n_point); % 计算样本数量
- 将
x转换为列向量，确保处理时的统一性。 n_sample计算可以分割出的完整样本数
3. 核心处理：不同时频分析方法
（1）Mel 频谱图 (mel)
[S, cF, t] = melSpectrogram(xx, fs, 'WindowLength', WL, 'OverlapLength', OL, 'FFTLength', FL);
melSpectrogram 计算 Mel 频谱，并转换为 dB 形式后绘制图像。
（2）短时傅里叶变换 (stft)
[S, F, T, ~] = spectrogram(xx, window, overlap, nfft, fs);
imagesc(T, F, 20*log10(abs(S)));
spectrogram 计算 STFT，得到时间-频率图，并转换为 dB 形式。
（3）连续小波变换 (cwt)
coefs = cwt(xx, scal, wavename);
f = scal2frq(scal, wavename, 1/fs);
imagesc(t, f, abs(coefs));
cwt 计算小波变换，并转换为频率坐标。
（4）S 变换 (ST)
[st_matrix, st_t, st_freq] = st(xx);
imagesc(st_t, st_freq, abs(st_matrix));
st 计算 S 变换（时间-频率分解）。
（5）希尔伯特-黄变换 (HHT)
imf = emd(xx);
hht(imf, Fs);
emd(xx)计算经验模态分解 (EMD)。hht(imf, Fs)绘制 Hilbert-Huang 变换 (HHT)。
（6）格拉姆角场 (GASF)
GM = GAF(xx, '+', 0);
imagesc(GM);
GAF 计算 Gramian Angular Field (GAF)。
（7）递归图 (RP)
rp_mat = RecurrencePlot(xx, m, tau, eta, 0);
imshow(rp_mat, []);
计算递归图 (Recurrence Plot, RP)。
（8）同步压缩变换 (SST)
Ts = SST(xx);
imagesc(time, fre, abs(Ts));
SST 计算同步压缩变换 (Synchrosqueezing Transform)。
4. 结果保存
saveas(gcf, fullfile(file_path, [method '_' num2str(i) '.jpg']))
使用 saveas 将生成的图像保存为 jpg 格式。
axis_control.m
function axis_control(flag)
% 用以控制坐标参数
image_size = 400; % 图像大小
if flag==1
ylabel('Frequency (Hz)')
xlabel('x')
else
% 无边，无框，无坐标轴
box off; axis off % 无框
set(gcf,'Position',[600 100 image_size image_size]) % 图像大小
set(gca,'xticklabel',[],'yticklabel',[]); % 无坐标轴
set(gca,'position',[0 0 1 1]) % 无边
end
endAlgorithms 文件夹
Wsst2_new.m
function [WSST2, fs, as, omega, omega2, tau, phipp, norm2psi] = Wsst2_new(s,gamma,mywav,nv)
% function Wsst : computes the synchrosqueezing transform of signal s.
% New formulae Cphi for reconstruction
%
% Inputs
% s : input signal, power of 2
% dt : sample period
% gamma : threshold
% mywav : string coding the mother wavelet
% nv : number of coefficient per octave
% doplot : 0 or 1 for extern plotting
%
% Outputs :
% WSST2: the second order % 基于小波二阶同步压缩变换
% fs : frequency vector
% as : scales vector
% nargin/out
if nargin<4
nv = 32; % nbre coef s par octave
gamma = 0.01; % Seuil pour wavelet thresholding
mywav = 'cmor2-1';
end
% Computing scales
s = s(:);
n = length(s);
noct = log2(n)-1;
na = noct*nv+1;
as = (2^(-1/nv)) .^ (0:1:na-1);
%noct = log2(n)-1;
%na = noct*nv;
%as = (2^(1/nv) .^ (1:1:na));
% Padding signal (symmetric)
[N, sx, n1] = mypad(s);
%% Wavelet transform and reassignment operators
WT = zeros(na, N); % CWT
omega = zeros(na, N); % Vertical operator
tau = zeros(na, N); % Horizontal one
phipp = zeros(na, N); % FM
omega2 = zeros(na, N); % Second-order IF
Dtau = zeros(na, N);
sx = sx(:).';
xh = fft(sx);
%xi = (0:N-1)/N;
xi = (0:N-1)/2; % to change if signal not a power of 2
% Filter definition
if strncmp(mywav,'gmor',4)
[v1 v2] = regexp(mywav,'[0-9.]*-[0-9.]');
beta = str2num(mywav(v1:v2-2));
gam = str2num(mywav(v2:end));
filt = @(a) gmor(beta,gam,a*xi);
func = @(x) 2*(exp(1)*gam/beta)^(beta/gam) * x.^beta .* exp(-x.^gam);
elseif strncmp(mywav,'cmor',4)
[v1 v2] = regexp(mywav,'[0-9.]*-[0-9.]');
Fb = str2num(mywav(v1:v2-2));
Fc= str2num(mywav(v2:end));
filt = @(a) cmor(Fb,Fc,a*xi);
func = @(x) sqrt(Fb) * exp(-Fb^2*pi*(x-Fc).^2);
elseif strncmp(mywav,'bump',4)
[v1 v2] = regexp(mywav,'[0-9.]*-[0-9.]');
mu = str2num(mywav(v2:end));
sigma = str2num(mywav(v1:v2-2));
filt = @(a) bump(mu,sigma,a*xi);
func = @(x) exp(1-1./(1-((x-mu)./sigma).^2));
end
% for each octave
norm2psi = zeros(1,length(as));
for ai = 1:length(as)
a = as(ai);
[psih, dfilt,~] = filt(a);
norm2psi(ai) = norm(psih);
% temporary values
Wtmp = ifft(conj(psih).* xh);
Wnu = ifft(conj(a*xi.*psih) .* xh);
Wnunu = ifft(conj((a*xi).^2.*psih) .* xh);
Wp = ifft(conj(dfilt).*xh);
Wnup = ifft(conj(a*xi.*dfilt).*xh);
% store results
WT(ai,:) = Wtmp;
omega(ai,:) = 1/a*Wnu./Wtmp;
tau(ai, :) = a*Wp./Wtmp/2/1i/pi;
phipp(ai,:) = 2*1i*pi/a^2*(Wnunu.*Wtmp-Wnu.^2)./(Wtmp.^2+Wnup.*Wtmp-Wp.*Wnu);
Dtau(ai,:) = (Wtmp.^2+Wnup.*Wtmp-Wp.*Wnu)./(Wtmp.^2);
end
%phipp(abs(real(tau)*n)<1)=0;
omega2 = real(omega - phipp.*tau);
omega = real(omega);
tau= real(tau);
WT = WT(:, n1+1:n1+n);
omega = omega(:, n1+1:n1+n);
tau = tau(:, n1+1:n1+n);
phipp = phipp(:, n1+1:n1+n);
omega2 = omega2(:, n1+1:n1+n);
Dtau = Dtau(:, n1+1:n1+n);
% figure();
% imagesc((0:n-1)/n,1./as,abs(Dtau)); set(gca,'YDir','normal');
%
% pause
%%
omega(omega<0) = NaN;
omega(abs(WT)<gamma) = NaN;
tau(abs(WT)<gamma) = NaN;
omega2(omega2<0) = NaN;
omega2(abs(WT)<gamma) = NaN;
phipp(abs(WT)<gamma) = NaN;
%% Synchrosqueezing and reassignment
fs = 1./as;
WSST = zeros(size(WT));
WSST2 = zeros(size(WT));
for b=1:n
for ai=1:length(as)
if (~isnan(omega(ai,b)) && abs(WT(ai,b))> 2*gamma*sqrt(2)*norm2psi(ai)/sqrt(2*n)) %module
%if (abs(Dtau(ai,b)) > 0)
k = round(1+nv*log2(omega2(ai,b)));
if k>0 && k<=na
WSST2(k,b) = WSST2(k, b) + WT(ai, b); % WSST2
end
% else
% k = round(1+nv*log2(omega(ai,b)));
% if k>0 && k<=na
% WSST2(k,b) = WSST2(k, b) + WT(ai, b); % WSST
% end
% end
k = round(1+nv*log2(omega(ai,b)));
if k>0 && k<=na
WSST(k,b) = WSST(k, b) + WT(ai, b); % WSST
end
end
end % for ai...
end % for b
WSST = WSST/nv*log(2);
WSST2 = WSST2/nv*log(2);
end
function [psih, dpsih, ddpsih] = cmor(Fb,Fc,xi)
psih = sqrt(Fb) * exp(-Fb^2*pi*(xi-Fc).^2);
dpsih = -2*Fb^2*pi*(xi-Fc).*psih;
%ddpsih = (-2*Fb^2*pi + (-2*Fb^2*pi*(xi-Fc)).^2).*psih;
ddpsih = 2*pi*Fb^2*(2*pi*Fb^2*Fc^2-4*pi*Fb^2*Fc.*xi+2*pi*Fb^2*xi.^2-1).*psih;
end
function [N, x, n1] = mypad(s)
% mypad : symmetric padding of vector s.
n = length(s);
% Computing the power of 2
N = 2^(1+round(log2(n+eps)));
n1 = floor((N-n)/2);
% Padding
padtype = 'symmetric';
try
sleft = padarray(s, n1, padtype, 'pre');
sright = padarray(s, n1, padtype, 'post');
catch
warning('Image processing Toolbox is not working. Using constant pading');
sleft = flipud(s(2:n1+1));
sright = flipud(s(end-n1:end-1));
end
x = [sleft(1:n1); s; sright(end-n1+1:end)];
endwset.m
function [sst,f] = wset(x,varargin)
%Wavelet Synchroextracted Transform
% SET = wset(X) returns the wavelet synchroextracted transform for the 1-D
% real-valued signal, X. X must have at least 4 samples. The
% synchroextracted transform uses 32 voices per octave and the number of
% octaves is floor(log2(numel(X)))-1. The transform uses the analytic
% Morlet wavelet by default. SET is a Na-by-N matrix where Na is the
% number of scales, 32*(floor(log2(numel(X)))-1), and N is the number of
% samples in X.
%
% [SET,F] = wset(X) returns a vector of frequencies, F, in cycles/sample
% corresponding to the rows of SET.
%
% [...] = wset(X,Fs) specifies the sampling frequency, Fs, in hertz as a
% positive scalar. If you specify the sampling frequency, WSET returns
% the frequencies in hertz.
%
% [...] = wset(X,Ts) uses the positive <a href="matlab:help duration">duration</a>, Ts, to compute the
% scale-to-frequency conversion, F. Ts is the time between samples of X.
% F has units of cycles/unit time where the unit of time is the same time
% unit as the duration.
%
% [...] = wset(...,WAV) uses the analytic wavelet specified by WAV to
% compute the synchroextracted transform. Valid choices for WAV are
% 'amor' and 'bump' for the analytic Morlet and bump wavelet. If
% unspecified, WAV defaults to 'amor'.
%
% [...] = wset(...,'VoicesPerOctave',NV) specifies the number of voices
% per octave as a positive even integer between 10 and 48. The number of
% scales is the product of the number of voices per octave and the number
% of octaves. If unspecified, NV defaults to 32 voices per octave. You
% can specify the 'VoicesPerOctave' name-value pair anywhere in the input
% argument list after the signal X.
%
% [...] = wset(...,'ExtendSignal',EXTENDFLAG) specifies whether to
% symmetrically extend the signal by reflection to mitigate boundary
% effects. EXTENDFLAG can be one of the following options [ true |
% {false}]. If unspecified, EXTENDSIGNAL defaults to false. You can
% specify the 'ExtendSignal' name-value pair anywhere in the input
% argument list after the signal X.
%
% wset(...) with no output arguments plots the wavelet synchroextracted
% transform as a function of time and frequency. If you do not specify a
% sampling frequency or interval, the synchroextracted transform is
% plotted in cycles/sample. If you supply a sampling frequency, Fs, the
% synchroextracted transform is plotted in hertz. If you supply a <a href="matlab:help duration">duration</a>
% as a sampling interval, the synchroextracted transform is plotted
% in cycles/unit time where the time unit is the same as the duration.
%
% % Example 1:
% % Obtain the wavelet synchroextracted transform of a quadratic chirp.
% % The chirp is sampled at 1000 Hz.
% load quadchirp;
% [set,f] = wset(quadchirp,1000);
% hp = pcolor(tquad,f,abs(set));
% hp.EdgeColor = 'none';
% title('Wavelet Synchroextracted Transform');
% xlabel('Time'); ylabel('Hz');
%
% % Example 2:
% % Obtain the wavelet synchroextracted transform of the sunspot
% % data. Specify the sampling interval to be 1 for one sample per
% % year.
% load sunspot.dat;
% wset(sunspot(:,2),years(1))
%
% See also iwsst, wsstridge, duration
narginchk(1,8);
nbSamp = numel(x);
x = x(:)';
validateattributes(x,{'double'},{'row','finite','real'},'wsst','X');
if numel(x)<4
error(message('Wavelet:synchroextracte:NumInputSamples'));
end
params = parseinputs(nbSamp,varargin{:});
nv = params.nv;
noct = params.noct;
% Create scale vector
na = noct*params.nv;
% If sampling frequency is specified, dt = 1/fs
if (isempty(params.fs) && isempty(params.Ts))
% The default is 1 for normalized frequency
dt = params.dt;
Units = '';
elseif (~isempty(params.fs) && isempty(params.Ts))
% Accept the sampling frequency in hertz
fs = params.fs;
dt = 1/fs;
Units = '';
elseif (isempty(params.fs) && ~isempty(params.Ts))
% Get the dt and Units from the duration object
[dt,Units] = getDurationandUnits(params.Ts);
end
a0 = 2^(1/nv);
scales = a0.^(1:na);
NbSc = numel(scales);
% Construct time series to analyze, pad if necessary
meanSIG = mean(x);
x = x - meanSIG;
NumExten = 0;
if params.pad
%Pad the time series symmetrically
np2 = nextpow2(nbSamp);
NumExten = 2^np2-nbSamp;
x = wextend('1d','symw',x,NumExten,'b');
end
%Record data length plus any extension
N = numel(x);
%Create frequency vector for CWT computation
omega = (1:fix(N/2));
omega = omega.*((2.*pi)/N);
omega = [0., omega, -omega(fix((N-1)/2):-1:1)];
% Compute FFT of the (padded) time series
xdft = fft(x);
[psift,dpsift] = sstwaveft(params.WAV,omega,scales,params.wavparam);
%Obtain CWT coefficients and derivative
cwtcfs = ifft(repmat(xdft,NbSc,1).*psift,[],2);
dcwtcfs = ifft(repmat(xdft,NbSc,1).*dpsift,[],2);
%Remove padding if any
cwtcfs = cwtcfs(:,NumExten+1:end-NumExten);
dcwtcfs = dcwtcfs(:,NumExten+1:end-NumExten);
%Compute the phase transform
phasetf = imag(dcwtcfs./cwtcfs)./(2*pi);
% Threshold for synchroextracting
phasetf(abs(phasetf)<params.thr) = NaN;
%Wavelet Synchroextracting Transform (WSET)
phasetf=round(phasetf*N);
phasetf((phasetf)<1) = 1;
phasetf((phasetf)>N/2) = 1;
phasetf(isnan(phasetf)) = 1;
[MM,NN]=size(phasetf);
phasescale=zeros(MM,NN);
fre2scale=zeros(1,round(NN/2));
for nn=1:round(NN/2)
[~,fre2scale(nn)]=max(psift(:,nn));
end
for mm=1:MM
for nn=1:NN
phasescale(mm,nn)= fre2scale(phasetf(mm,nn));
end
end
cwtcfs2=zeros(size(cwtcfs));
for mm=1:MM
for nn=1:NN
if phasescale(mm,nn)==mm
cwtcfs2(mm,nn)= cwtcfs(mm,nn);
end
end
end
cwtcfs=cwtcfs2;
phasetf=phasetf/N;
%
% Create frequency vector for output
log2Nyquist = log2(1/(2*dt));
log2Fund = log2(1/(nbSamp*dt));
freq = 2.^linspace(log2Fund,log2Nyquist,na);
Tx = 1/nv*sstalgo(cwtcfs,phasetf,params.thr);
if (nargout == 0)
plotsst(Tx,freq,dt,params.engunitflag,params.normalizedfreq,Units);
else
sst = Tx;
f = freq;
end
%-------------------------------------------------------------------------
function [wft,dwft] = sstwaveft(WAV,omega,scales,wavparam)
% Admissible wavelets are:
% - MORLET wavelet (A) - 'morl':
% PSI_HAT(s) = exp(-(s-s0).^2/2) * (s>0)
% Parameter: s0, default s0 = 6.
% - Bump wavelet: 'bump':
% PSI_HAT(s) = exp(1-(1/((s-mu)^2./sigma^2))).*(abs((s-mu)/sigma)<1)
% Parameters: mu,sigma.
% default: mu=5, sigma = 0.6.
% Normalized to have unit magnitude at the peak frequency of the wavelet
NbSc = numel(scales);
NbFrq = numel(omega);
wft = zeros(NbSc,NbFrq);
switch WAV
case 'amor'
cf = wavparam;
for jj = 1:NbSc
expnt = -(scales(jj).*omega - cf).^2/2.*(omega > 0);
wft(jj,:) = exp(expnt).*(omega > 0);
end
case 'bump'
mu = wavparam(1);
sigma = wavparam(2);
for jj = 1:NbSc
w = (scales(jj)*omega-mu)./sigma;
expnt = -1./(1-w.^2);
daughter = exp(1)*exp(expnt).*(abs(w)<1-eps(1));
daughter(isnan(daughter)) = 0;
wft(jj,:) = daughter;
end
end
%Compute derivative
omegaMatrix = repmat(omega,NbSc,1);
dwft = 1j*omegaMatrix.*wft;
%-------------------------------------------------------------------------
function plotsst(Tx,F,dt,engunitflag,isfreqnormalized,Units)
if ~isempty(Units)
freqUnits = Units(1:end-1);
end
t = 0:dt:(size(Tx,2)*dt)-dt;
if engunitflag && isfreqnormalized
frequnitstrs = wgetfrequnitstrs;
freqlbl = frequnitstrs{1};
xlbl = 'Samples';
elseif engunitflag && ~isfreqnormalized
[F,~,uf] = engunits(F,'unicode');
freqlbl = wgetfreqlbl([uf 'Hz']);
[t,~,ut] = engunits(t,'unicode','time');
xlbl = [getString(message('Wavelet:getfrequnitstrs:Time')) ' (' ut ')'];
else
freqlbl = getString(message('Wavelet:synchrosqueezed:FreqLabel'));
freqlbl = ...
[freqlbl '/' freqUnits ')'];
xlbl = getString(message('Wavelet:synchrosqueezed:Time'));
xlbl = [xlbl ' (' Units ')'];
end
h = pcolor(t,F,abs(Tx));
h.EdgeColor = 'none';
shading interp;
ylabel(freqlbl); xlabel(xlbl);
title(getString(message('Wavelet:synchrosqueezed:SynchrosqueezedTitle')));
%-------------------------------------------------------------------------
function params = parseinputs(nbSamp,varargin)
% Set defaults.
params.fs = [];
params.dt = 1;
params.Ts = [];
params.sampinterval = false;
params.engunitflag = true;
params.WAV = 'amor';
params.wavparam = 6;
params.thr = 1e-8;
params.nv = 32;
params.noct = floor(log2(nbSamp))-1;
params.pad = false;
params.normalizedfreq = true;
% Error out if there are any calendar duration objects
tfcalendarDuration = cellfun(@iscalendarduration,varargin);
if any(tfcalendarDuration)
error(message('Wavelet:FunctionInput:CalendarDurationSupport'));
end
tfsampinterval = cellfun(@isduration,varargin);
if (any(tfsampinterval) && nnz(tfsampinterval) == 1)
params.sampinterval = true;
params.Ts = varargin{tfsampinterval>0};
if (numel(params.Ts) ~= 1 ) || params.Ts <= 0 || isempty(params.Ts)
error(message('Wavelet:FunctionInput:PositiveScalarDuration'));
end
params.engunitflag = false;
params.normalizedfreq = false;
varargin(tfsampinterval) = [];
end
%Look for Name-Value pairs
numvoices = find(strncmpi('voicesperoctave',varargin,1));
if any(numvoices)
params.nv = varargin{numvoices+1};
%validate the value is logical
validateattributes(params.nv,{'numeric'},{'positive','scalar',...
'even','>=',10,'<=',48},'wsst','VoicesPerOctave');
varargin(numvoices:numvoices+1) = [];
if isempty(varargin)
return;
end
end
extendsignal = find(strncmpi('extendsignal',varargin,1));
if any(extendsignal)
params.pad = varargin{extendsignal+1};
if ~isequal(params.pad,logical(params.pad))
error(message('Wavelet:FunctionInput:Logical'));
end
varargin(extendsignal:extendsignal+1) = [];
if isempty(varargin)
return;
end
end
% Only scalar left must be sampling frequency or sampling interval
% Only scalar left must be sampling frequency
tfsampfreq = cellfun(@(x) (isscalar(x) && isnumeric(x)),varargin);
if (any(tfsampfreq) && (nnz(tfsampfreq) == 1) && ~params.sampinterval)
params.fs = varargin{tfsampfreq};
validateattributes(params.fs,{'numeric'},{'positive'},'wsst','Fs');
params.normalizedfreq = false;
params.engunits = true;
elseif any(tfsampfreq) && params.sampinterval
error(message('Wavelet:FunctionInput:SamplingIntervalOrDuration'));
elseif nnz(tfsampfreq)>1
error(message('Wavelet:FunctionInput:Invalid_ScalNum'));
end
%Only char variable left must be wavelet
tfwav = cellfun(@ischar,varargin);
if (nnz(tfwav) == 1)
params.WAV = varargin{tfwav>0};
params.WAV = validatestring(params.WAV,{'bump','amor'},'wsst','WAV');
elseif nnz(tfwav)>1
error(message('Wavelet:FunctionInput:InvalidChar'));
end
if strncmpi(params.WAV,'bump',1)
params.wavparam = [5 1];
end
%------------------------------------------------------------------------
function Tx = sstalgo(cwtcfs,phasetf,gamma)
M = size(cwtcfs,1);
N = size(cwtcfs,2);
log2Fund = log2(1/N);
log2Nyquist = log2(1/2);
iRow = real(1 + floor(M/(log2Nyquist-log2Fund)*(log2(phasetf)-log2Fund)));
idxphasetf = find(iRow>0 & iRow<=M & ~isnan(iRow));
idxcwtcfs = find(abs(cwtcfs)>gamma);
idx = intersect(idxphasetf,idxcwtcfs);
iCol = repmat(1:N,M,1);
Tx = accumarray([iRow(idx) iCol(idx)],cwtcfs(idx),size(cwtcfs));
wmsst.m
function [sst,f] = wmsst(x,num,varargin)
%Wavelet Multisynchrosqueezed Transform
% MSST = wmsst(X) returns the wavelet multisynchrosqueezed transform for the 1-D
% real-valued signal, X. X must have at least 4 samples. The
% multisynchrosqueezed transform uses 32 voices per octave and the number of
% octaves is floor(log2(numel(X)))-1. The transform uses the analytic
% Morlet wavelet by default. MSST is a Na-by-N matrix where Na is the
% number of scales, 32*(floor(log2(numel(X)))-1), and N is the number of
% samples in X.
%
% [SST,F] = wmsst(X,Num) returns a vector of frequencies, F, in cycles/sample
% corresponding to the rows of MSST.
%
% [...] = wmsst(X,Num,Fs) specifies the sampling frequency, Fs, in hertz as a
% positive scalar and Num is the iterative number. If you specify the sampling frequency, WMSST returns
% the frequencies in hertz.
%
% [...] = wmsst(X,Num,Ts) uses the positive <a href="matlab:help duration">duration</a>, Ts, to compute the
% scale-to-frequency conversion, F. Ts is the time between samples of X.
% F has units of cycles/unit time where the unit of time is the same time
% unit as the duration.Num is the iterative number.
%
% [...] = wmsst(...,WAV) uses the analytic wavelet specified by WAV to
% compute the multisynchrosqueezed transform. Valid choices for WAV are
% 'amor' and 'bump' for the analytic Morlet and bump wavelet. If
% unspecified, WAV defaults to 'amor'.
%
% [...] = wmsst(...,'VoicesPerOctave',NV) specifies the number of voices
% per octave as a positive even integer between 10 and 48. The number of
% scales is the product of the number of voices per octave and the number
% of octaves. If unspecified, NV defaults to 32 voices per octave. You
% can specify the 'VoicesPerOctave' name-value pair anywhere in the input
% argument list after the signal X.
%
% [...] = wmsst(...,'ExtendSignal',EXTENDFLAG) specifies whether to
% symmetrically extend the signal by reflection to mitigate boundary
% effects. EXTENDFLAG can be one of the following options [ true |
% {false}]. If unspecified, EXTENDSIGNAL defaults to false. You can
% specify the 'ExtendSignal' name-value pair anywhere in the input
% argument list after the signal X.
%
% wmsst(...) with no output arguments plots the wavelet multisynchrosqueezed
% transform as a function of time and frequency. If you do not specify a
% sampling frequency or interval, the multisynchrosqueezed transform is
% plotted in cycles/sample. If you supply a sampling frequency, Fs, the
% multisynchrosqueezed transform is plotted in hertz. If you supply a <a href="matlab:help duration">duration</a>
% as a sampling interval, the multisynchrosqueezed transform is plotted
% in cycles/unit time where the time unit is the same as the duration.
%
% % Example 1:
% % Obtain the wavelet multisynchrosqueezed transform of a quadratic chirp.
% % The chirp is sampled at 1000 Hz.
% load quadchirp;
% [msst,f] = wmsst(quadchirp,10,1000);
% hp = pcolor(tquad,f,abs(msst));
% hp.EdgeColor = 'none';
% title('Wavelet Multisynchrosqueezed Transform');
% xlabel('Time'); ylabel('Hz');
%
% % Example 2:
% % Obtain the wavelet multisynchrosqueezed transform of the sunspot
% % data. Specify the sampling interval to be 1 for one sample per
% % year.
% load sunspot.dat;
% wmsst(sunspot(:,2),10,years(1))
%
% See also iwsst, wsstridge, duration
narginchk(1,9);
nbSamp = numel(x);
x = x(:)';
validateattributes(x,{'double'},{'row','finite','real'},'wsst','X');
if numel(x)<4
error(message('Wavelet:synchrosqueezed:NumInputSamples'));
end
params = parseinputs(nbSamp,varargin{:});
nv = params.nv;
noct = params.noct;
% Create scale vector
na = noct*params.nv;
% If sampling frequency is specified, dt = 1/fs
if (isempty(params.fs) && isempty(params.Ts))
% The default is 1 for normalized frequency
dt = params.dt;
Units = '';
elseif (~isempty(params.fs) && isempty(params.Ts))
% Accept the sampling frequency in hertz
fs = params.fs;
dt = 1/fs;
Units = '';
elseif (isempty(params.fs) && ~isempty(params.Ts))
% Get the dt and Units from the duration object
[dt,Units] = getDurationandUnits(params.Ts);
end
a0 = 2^(1/nv);
scales = a0.^(1:na);
NbSc = numel(scales);
% Construct time series to analyze, pad if necessary
meanSIG = mean(x);
x = x - meanSIG;
NumExten = 0;
if params.pad
%Pad the time series symmetrically
np2 = nextpow2(nbSamp);
NumExten = 2^np2-nbSamp;
x = wextend('1d','symw',x,NumExten,'b');
end
%Record data length plus any extension
N = numel(x);
%Create frequency vector for CWT computation
omega = (1:fix(N/2));
omega = omega.*((2.*pi)/N);
omega = [0., omega, -omega(fix((N-1)/2):-1:1)];
% Compute FFT of the (padded) time series
xdft = fft(x);
[psift,dpsift] = sstwaveft(params.WAV,omega,scales,params.wavparam);
%Obtain CWT coefficients and derivative
cwtcfs = ifft(repmat(xdft,NbSc,1).*psift,[],2);
dcwtcfs = ifft(repmat(xdft,NbSc,1).*dpsift,[],2);
%Remove padding if any
cwtcfs = cwtcfs(:,NumExten+1:end-NumExten);
dcwtcfs = dcwtcfs(:,NumExten+1:end-NumExten);
%Compute the phase transform
phasetf = imag(dcwtcfs./cwtcfs)./(2*pi);
% Threshold for synchrosqueezing
phasetf(abs(phasetf)<params.thr) = NaN;
%Wavelet Multisynchrosqueezing Transform (WMSST)
if num>1
phasetf=round(phasetf*N);
phasetf((phasetf)<1) = 1;
phasetf((phasetf)>N/2) = 1;
phasetf(isnan(phasetf)) = 1;
[MM,NN]=size(phasetf);
phasescale=zeros(MM,NN);
phasescale2=zeros(MM,NN);
phasetf2=zeros(MM,NN);
fre2scale=zeros(1,round(NN/2));
for nn=1:round(NN/2)
[~,fre2scale(nn)]=max(psift(:,nn));
end
scale2fre=zeros(1,MM);
for mm=1:MM
[~,scale2fre(mm)]=max(psift(mm,:));
end
for mm=1:MM
for nn=1:NN
phasescale(mm,nn)= fre2scale(phasetf(mm,nn));
end
end
for kk=1:num-1
for nn=1:NN
for mm=1:MM
k = phasescale(mm,nn);
if k>=1 && k<=NN
phasescale2(mm,nn)=phasescale(k,nn);
end
end
end
phasescale=phasescale2;
end
for mm=1:MM
for nn=1:NN
phasetf2(mm,nn)= scale2fre(phasescale2(mm,nn));
end
end
phasetf=phasetf2/N;
end
% Create frequency vector for output
log2Nyquist = log2(1/(2*dt));
log2Fund = log2(1/(nbSamp*dt));
freq = 2.^linspace(log2Fund,log2Nyquist,na);
Tx = 1/nv*sstalgo(cwtcfs,phasetf,params.thr);
if (nargout == 0)
plotsst(Tx,freq,dt,params.engunitflag,params.normalizedfreq,Units);
else
sst = Tx;
f = freq;
end
%-------------------------------------------------------------------------
function [wft,dwft] = sstwaveft(WAV,omega,scales,wavparam)
% Admissible wavelets are:
% - MORLET wavelet (A) - 'morl':
% PSI_HAT(s) = exp(-(s-s0).^2/2) * (s>0)
% Parameter: s0, default s0 = 6.
% - Bump wavelet: 'bump':
% PSI_HAT(s) = exp(1-(1/((s-mu)^2./sigma^2))).*(abs((s-mu)/sigma)<1)
% Parameters: mu,sigma.
% default: mu=5, sigma = 0.6.
% Normalized to have unit magnitude at the peak frequency of the wavelet
NbSc = numel(scales);
NbFrq = numel(omega);
wft = zeros(NbSc,NbFrq);
switch WAV
case 'amor'
cf = wavparam;
for jj = 1:NbSc
expnt = -(scales(jj).*omega - cf).^2/2.*(omega > 0);
wft(jj,:) = exp(expnt).*(omega > 0);
end
case 'bump'
mu = wavparam(1);
sigma = wavparam(2);
for jj = 1:NbSc
w = (scales(jj)*omega-mu)./sigma;
expnt = -1./(1-w.^2);
daughter = exp(1)*exp(expnt).*(abs(w)<1-eps(1));
daughter(isnan(daughter)) = 0;
wft(jj,:) = daughter;
end
end
%Compute derivative
omegaMatrix = repmat(omega,NbSc,1);
dwft = 1j*omegaMatrix.*wft;
%-------------------------------------------------------------------------
function plotsst(Tx,F,dt,engunitflag,isfreqnormalized,Units)
if ~isempty(Units)
freqUnits = Units(1:end-1);
end
t = 0:dt:(size(Tx,2)*dt)-dt;
if engunitflag && isfreqnormalized
frequnitstrs = wgetfrequnitstrs;
freqlbl = frequnitstrs{1};
xlbl = 'Samples';
elseif engunitflag && ~isfreqnormalized
[F,~,uf] = engunits(F,'unicode');
freqlbl = wgetfreqlbl([uf 'Hz']);
[t,~,ut] = engunits(t,'unicode','time');
xlbl = [getString(message('Wavelet:getfrequnitstrs:Time')) ' (' ut ')'];
else
freqlbl = getString(message('Wavelet:synchrosqueezed:FreqLabel'));
freqlbl = ...
[freqlbl '/' freqUnits ')'];
xlbl = getString(message('Wavelet:synchrosqueezed:Time'));
xlbl = [xlbl ' (' Units ')'];
end
h = pcolor(t,F,abs(Tx));
h.EdgeColor = 'none';
shading interp;
ylabel(freqlbl); xlabel(xlbl);
title(getString(message('Wavelet:synchrosqueezed:SynchrosqueezedTitle')));
%-------------------------------------------------------------------------
function params = parseinputs(nbSamp,varargin)
% Set defaults.
params.fs = [];
params.dt = 1;
params.Ts = [];
params.sampinterval = false;
params.engunitflag = true;
params.WAV = 'amor';
params.wavparam = 6;
params.thr = 1e-8;
params.nv = 32;
params.noct = floor(log2(nbSamp))-1;
params.pad = false;
params.normalizedfreq = true;
% Error out if there are any calendar duration objects
tfcalendarDuration = cellfun(@iscalendarduration,varargin);
if any(tfcalendarDuration)
error(message('Wavelet:FunctionInput:CalendarDurationSupport'));
end
tfsampinterval = cellfun(@isduration,varargin);
if (any(tfsampinterval) && nnz(tfsampinterval) == 1)
params.sampinterval = true;
params.Ts = varargin{tfsampinterval>0};
if (numel(params.Ts) ~= 1 ) || params.Ts <= 0 || isempty(params.Ts)
error(message('Wavelet:FunctionInput:PositiveScalarDuration'));
end
params.engunitflag = false;
params.normalizedfreq = false;
varargin(tfsampinterval) = [];
end
%Look for Name-Value pairs
numvoices = find(strncmpi('voicesperoctave',varargin,1));
if any(numvoices)
params.nv = varargin{numvoices+1};
%validate the value is logical
validateattributes(params.nv,{'numeric'},{'positive','scalar',...
'even','>=',10,'<=',48},'wsst','VoicesPerOctave');
varargin(numvoices:numvoices+1) = [];
if isempty(varargin)
return;
end
end
extendsignal = find(strncmpi('extendsignal',varargin,1));
if any(extendsignal)
params.pad = varargin{extendsignal+1};
if ~isequal(params.pad,logical(params.pad))
error(message('Wavelet:FunctionInput:Logical'));
end
varargin(extendsignal:extendsignal+1) = [];
if isempty(varargin)
return;
end
end
% Only scalar left must be sampling frequency or sampling interval
% Only scalar left must be sampling frequency
tfsampfreq = cellfun(@(x) (isscalar(x) && isnumeric(x)),varargin);
if (any(tfsampfreq) && (nnz(tfsampfreq) == 1) && ~params.sampinterval)
params.fs = varargin{tfsampfreq};
validateattributes(params.fs,{'numeric'},{'positive'},'wsst','Fs');
params.normalizedfreq = false;
params.engunits = true;
elseif any(tfsampfreq) && params.sampinterval
error(message('Wavelet:FunctionInput:SamplingIntervalOrDuration'));
elseif nnz(tfsampfreq)>1
error(message('Wavelet:FunctionInput:Invalid_ScalNum'));
end
%Only char variable left must be wavelet
tfwav = cellfun(@ischar,varargin);
if (nnz(tfwav) == 1)
params.WAV = varargin{tfwav>0};
params.WAV = validatestring(params.WAV,{'bump','amor'},'wsst','WAV');
elseif nnz(tfwav)>1
error(message('Wavelet:FunctionInput:InvalidChar'));
end
if strncmpi(params.WAV,'bump',1)
params.wavparam = [5 1];
end
%------------------------------------------------------------------------
function Tx = sstalgo(cwtcfs,phasetf,gamma)
M = size(cwtcfs,1);
N = size(cwtcfs,2);
log2Fund = log2(1/N);
log2Nyquist = log2(1/2);
iRow = real(1 + floor(M/(log2Nyquist-log2Fund)*(log2(phasetf)-log2Fund)));
idxphasetf = find(iRow>0 & iRow<=M & ~isnan(iRow));
idxcwtcfs = find(abs(cwtcfs)>gamma);
idx = intersect(idxphasetf,idxcwtcfs);
iCol = repmat(1:N,M,1);
Tx = accumarray([iRow(idx) iCol(idx)],cwtcfs(idx),size(cwtcfs));
VSST2.m
function [VSST,omega,tau,omega2,g] = VSST2(s,aa,Nfft,gamma)
%% vertical second-order synchrosqueezing 垂直二阶同步压缩变换
%
% INPUTS:
% s: real or complex signal
% aa: the parameter a in the amgauss function of the Gaussian window
% Nfft: number of frequency bins
% gamma: threshold on the STFT for reassignment
% OUTPUTS:
% VSST : vertical second-order synchrosqueezing
% REFERENCES
% [1] Pham, D-H., Meignen, S. Second-order Synchrosqueezing Transform: The
% Wavelet Case and Comparisons
% nargin/out
if nargin<4
aa = 64; %
gamma = 0.001;
Nfft =1024;
end
s = s(:);
ft = 1:Nfft;
bt = 1:length(s);
nb = length(bt);
neta = length(ft);
prec = 10^(-6) ;
L = sqrt(Nfft/aa);
l = floor(sqrt(-Nfft*log(prec)/(aa*pi)))+1;
g = amgauss(2*l+1,l+1,L);
%figure(); plot(g)
% Window definition
n = (0:2*l)'-l;
t0 = n/Nfft;
t0 = t0(:);
a = aa*Nfft*pi;
gp = -2*a*t0.*g;
gpp = (-2*a+4*a^2*t0.^2).*g; % g''
% Initialization
STFT = zeros(neta,nb);
SST = zeros(neta,nb);
VSST = zeros(neta,nb);
omega = zeros(neta,nb);
tau = zeros(neta,nb);
omega2 = zeros(neta,nb);
phipp = zeros(neta,nb);
%% Computes STFT and reassignment operators
for b=1:nb
% STFT, window g
time_inst = -min([l,bt(b)-1]):min([l,nb-bt(b)]);
tmp(1:length(time_inst)) = s(bt(b)+time_inst).*g(l+1+time_inst);
A = fft(tmp(:),Nfft);
STFT(:,b) = A.*exp(-2/Nfft*pi*1i*(0:Nfft-1)'*time_inst(1))/Nfft*g(l+1);%renormalized so that it fits with recmodes
vg = A;
% STFT, window xg
tmp(1:length(time_inst)) = s(bt(b)+time_inst).*(time_inst)'/Nfft.*g(l+1+time_inst);
vxg = fft(tmp(:),Nfft);
% operator Lx (dtau)
tau(:,b) = vxg./vg;
% STFT, window gp
tmp(1:length(time_inst))= s(bt(b)+time_inst).*gp(l+1+time_inst);
vgp = fft(tmp(:),Nfft);
omega(:,b) = (ft-1)'- real(vgp/2/1i/pi./vg);
% STFT, window gpp
tmp(1:length(time_inst)) = s(bt(b)+time_inst).*gpp(l+1+time_inst);
vgpp = fft(tmp(:),Nfft);
%STFT, windox xgp
tmp(1:length(time_inst)) = s(bt(b)+time_inst).*(time_inst)'/Nfft.*gp(l+1+time_inst);
vxgp = fft(tmp(:),Nfft);
%computation of the two different omega
phipp(:,b) = 1/2/1i/pi*(vgpp.*vg-vgp.^2)./(vxg.*vgp-vxgp.*vg);
%new omega2
omega2(:,b) = omega(:,b) - real(phipp(:,b)).*real(tau(:,b))...
+ imag(phipp(:,b)).*imag(tau(:,b));
%omega2(:,b) = omega(:,b) - real(phipp(:,b)).*real(tau(:,b));
end
%% reassignment step
for b=1:nb
for eta=1:neta
if (abs(STFT(eta,b))> 2*sqrt(2)*gamma*norm(g)/length(s))
k = 1+round(omega(eta,b));
if (k >= 1) && (k <= neta)
% original reassignment
SST(k,b) = SST(k,b) + STFT(eta,b);
end
k = 1+round(omega2(eta,b));
if k>=1 && k<=neta
% second-order Vertical reassignment: VSST
VSST(k,b) = VSST(k,b) + STFT(eta,b);
end
end
end
end
end
function y = amgauss(N,t0,T);
%AMGAUSS Generate gaussian amplitude modulation.
% Y=AMGAUSS(N,T0,T) generates a gaussian amplitude modulation
% centered on a time T0, and with a spread proportional to T.
% This modulation is scaled such that Y(T0)=1
% and Y(T0+T/2) and Y(T0-T/2) are approximately equal to 0.5 .
%
% N : number of points.
% T0 : time center (default : N/2).
% T : time spreading (default : 2*sqrt(N)).
% Y : signal.
%
% Examples:
% z=amgauss(160); plot(z);
% z=amgauss(160,90,40); plot(z);
% z=amgauss(160,180,50); plot(z);
%
% See also AMEXPO1S, AMEXPO2S, AMRECT, AMTRIANG.
% F. Auger, July 1995.
% Copyright (c) 1996 by CNRS (France).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
if (nargin == 0),
error ( 'The number of parameters must be at least 1.' );
elseif (nargin == 1),
t0=N/2; T=2*sqrt(N);
elseif (nargin ==2),
T=2*sqrt(N);
end;
if (N<=0),
error('N must be greater or equal to 1.');
else
tmt0=(1:N)'-t0;
y = exp(-(tmt0/T).^2 * pi);
end
endTMSST.m
function [Ts, tfr2] = TMSST(x,hlength,num)
% Time-reassigned Multisynchrosqueezing Transform
% Input:
% x : Signal.
% hlength : Window length.
% num : Iteration number
% Output:
% Ts : TMSST result.
% tfr2 : STFT result.
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
%
% Written by YuGang.
[xrow,xcol] = size(x);
N=xrow;
if (xcol~=1),
error('X must be column vector');
end;
if (nargin < 2)
hlength=round(xrow/5);
num=1;
else if (nargin < 3)
num=1;
end
end
t=1:N;tcol=N;
hlength=hlength+1-rem(hlength,2);
ht = linspace(-0.5,0.5,hlength);ht=ht';
% Gaussian window
h = exp(-pi/0.32^2*ht.^2);
th=h.*ht;
[hrow,~]=size(h); Lh=(hrow-1)/2;
tfr1= zeros (N,tcol);
tfr3= zeros (N,tcol);
for icol=1:tcol,
ti= t(icol); tau=-min([round(N/2)-1,Lh,ti-1]):min([round(N/2)-1,Lh,xrow-ti]);
indices= rem(N+tau,N)+1;
rSig = x(ti+tau,1);
tfr1(indices,icol)=rSig.*conj(h(Lh+1+tau));
tfr3(indices,icol)=rSig.*conj(th(Lh+1+tau));
end;
tfr1=fft(tfr1);
tfr3=fft(tfr3);
tfr1=tfr1(1:round(N/2),:);
tfr3=tfr3(1:round(N/2),:);
omega2= zeros(round(N/2),tcol);
omega22= zeros(round(N/2),tcol);
for a=1:round(N/2)
omega2(a,:) = t+(hlength-1)*real(tfr3(a,t)./tfr1(a,t));
end
[neta,nb]=size(tfr1);
if num>1
for kk=1:num-1
for b=1:nb
for eta=1:neta
k2 = round(omega2(eta,b));
if k2>=1 && k2<=nb
omega22(eta,b)=omega2(eta,k2);
end
end
end
omega2=omega22;
end
else
omega22=omega2;
end
omega22=round(round(omega22*2)/2);
tfr2 = zeros(round(N/2),tcol);
for eta=1:round(N/2)%frequency
tfr2(eta,:)=tfr1(eta,:).*exp(-1j * 2 * pi*eta*(t/N));
end
Ts = zeros(round(N/2),tcol);
% Reassignment step
for b=1:N%time
for eta=1:round(N/2)%frequency
% if abs(tfr1(eta,b))>0.000001
k2 = omega22(eta,b);
if k2>=1 && k2<=N
Ts(eta,k2) = Ts(eta,k2) + (tfr2(eta,b));
end
% end
end
end
Ts=Ts/(xrow/2);
TET.m
function [Te] = TET(x,hlength)
% Transient-extracting transform
% Input
% x : Signal.
% hlength : Window length.
% Output
% Te : TET Representation.
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
%
% Written at 2017.3.17.
[xrow,xcol] = size(x);
N=xrow;
if (xcol~=1),
error('X must be column vector');
end;
if (nargin < 2),
hlength=round(xrow/8);
end;
t=1:N;
%ft = 1:round(N/2);
[~,tcol] = size(t);
hlength=hlength+1-rem(hlength,2);
ht = linspace(-0.5,0.5,hlength);ht=ht';
% Gaussian window
h = exp(-pi/0.32^2*ht.^2);
% Window tg(t)
th=h.*ht;
[hrow,hcol]=size(h); Lh=(hrow-1)/2;
tfr1= zeros (N,tcol);
tfr2= zeros (N,tcol);
for icol=1:tcol,
ti= t(icol); tau=-min([round(N/2)-1,Lh,ti-1]):min([round(N/2)-1,Lh,xrow-ti]);
indices= rem(N+tau,N)+1;
rSig = x(ti+tau,1);
tfr1(indices,icol)=rSig.*conj(h(Lh+1+tau));
tfr2(indices,icol)=rSig.*conj(th(Lh+1+tau));
end;
tfr1=fft(tfr1);
tfr2=fft(tfr2);
tfr1=tfr1(1:round(N/2),:);
tfr2=tfr2(1:round(N/2),:);
E=mean(abs(x));
omega= zeros(round(N/2),tcol);
for a=1:round(N/2)
omega(a,:) = t+(hlength-1)*real(tfr2(a,t)./tfr1(a,t));
end
GD=omega;
TEO=zeros(round(N/2),tcol);
for i=1:round(N/2)%frequency
for j=1:N%time
if abs(tfr1(i,j))>12*E%
if abs(omega(i,j)-j)<0.5%
TEO(i,j)=1;
end
end
end
end
tfr=tfr1/(xrow/2);%the amplitude of tfr result has been pre-rectified.
%tfr=tfr1/(sum(h)/2);%
Te=TEO.*tfr;
STET.m
function [Te2] = STET(x,hlength)
% Computes the second-order transient-extracting transform (STET) of the signal x.
% INPUT
% x : Signal needed to be column vector.
% hlength: The length of window function.
% OUTPUT
% Te2 : The STET result
[xrow,xcol] = size(x);
if (xcol~=1),
error('X must be column vector');
end;
N=xrow;
% hlength: the length of window function.
hlength=hlength+1-rem(hlength,2);
ht = linspace(-0.5,0.5,hlength);ht=ht';
% Gaussian window
h = exp(-0.5*ht.^2/0.02);
% derivative of window
dh = -ht .* h/0.02; % g'
[hrow,~]=size(h); Lh=(hrow-1)/2;
tfr1= zeros (round(N/2),N) ;% G(t,w);
tfr2= zeros (round(N/2),N) ;% Gg'(t,w);
tfr3= zeros (round(N/2),N) ;% Gsu(t,w);
t=1:N;
su=x.*t';% s*u;
for icol=1:N,
ti= t(icol); tau=-min([round(N/2)-1,Lh,ti-1]):min([round(N/2)-1,Lh,xrow-ti]);
indices= rem(N+tau,N)+1;
rSig = x(ti+tau,1);
suSig = su(ti+tau,1);
tfr1(indices,icol)=rSig.*h(Lh+1+tau);% G(t,w)*exp(iwt);
tfr2(indices,icol)=rSig.*dh(Lh+1+tau);% Gg'(t,w)*exp(iwt);
tfr3(indices,icol)=suSig.*h(Lh+1+tau);% Gsu(t,w)*exp(iwt);
end;
tfr1=fft(tfr1);tfr1=tfr1(1:round(N/2),:);
tfr2=fft(tfr2);tfr2=tfr2(1:round(N/2),:);
tfr3=fft(tfr3);tfr3=tfr3(1:round(N/2),:);
for eta=1:round(N/2)
tfr1(eta,:)=tfr1(eta,:).*exp(-1j * 2 * pi*(eta-1)*((t)/N));% G(t,w);
tfr2(eta,:)=tfr2(eta,:).*exp(-1j * 2 * pi*(eta-1)*((t)/N));% Gg'(t,w);
tfr3(eta,:)=tfr3(eta,:).*exp(-1j * 2 * pi*(eta-1)*((t)/N));% Gsu(t,w);
end
omega_t=real(tfr3./tfr1);
omega_w=real(-tfr2./tfr1);
omega_t_dt = zeros (round(N/2),N-1);
omega_w_dt = zeros (round(N/2),N-1);
omega_t_dw = zeros (round(N/2)-1,N);
omega_w_dw = zeros (round(N/2)-1,N);
for i=1:round(N/2)
omega_t_dt(i,:)=diff(omega_t(i,:));
omega_w_dt(i,:)=diff(omega_w(i,:));
end
omega_t_dt(:,end+1)=omega_t_dt(:,end);
omega_w_dt(:,end+1)=omega_w_dt(:,end);
for i=1:N
omega_t_dw(:,i)=diff(omega_t(:,i));
omega_w_dw(:,i)=diff(omega_w(:,i));
end
omega_t_dw(end+1,:)=omega_t_dw(end,:);
omega_w_dw(end+1,:)=omega_w_dw(end,:);
omega_t2=zeros(round(N/2),N);
tfr=zeros(round(N/2),N);
tfr1=tfr1/(xrow/2);
for i=1:round(N/2)%frequency variable
for j=1:N %time variable
if abs(omega_w_dt(i,j).*omega_t_dw(i,j))>1e-10;
omega_t2(i,j)=-omega_w_dw(i,j)/(omega_w_dt(i,j)*omega_t_dw(i,j))*(omega_t(i,j)-t(j)*omega_t_dt(i,j));
tfr(i,j)=tfr1(i,j)/(sqrt(-omega_w_dt(i,j)*omega_t_dw(i,j)/omega_w_dw(i,j)-1i*omega_w_dw(i,j)));
else
omega_t2(i,j)=omega_t(i,j);
tfr(i,j)=tfr1(i,j);
end
end
end
Te2=zeros(round(N/2),N);
omega_t2=round(omega_t2);
for i=1:round(N/2)
for j=1:N
if abs(omega_t2(i,j)-j)==0
Te2(i,j)=tfr(i,j);
end
end
end
endst.m
function [st,t,f] = st(timeseries,minfreq,maxfreq,samplingrate,freqsamplingrate)
% Returns the Stockwell Transform of the timeseries.
% Code by Robert Glenn Stockwell.
% DO NOT DISTRIBUTE
% BETA TEST ONLY
% Reference is "Localization of the Complex Spectrum: The S Transform"
% from IEEE Transactions on Signal Processing, vol. 44., number 4, April 1996, pages 998-1001.
%
%-------Inputs Needed------------------------------------------------
%
% *****All frequencies in (cycles/(time unit))!******
% "timeseries" - vector of data to be transformed
%-------Optional Inputs ------------------------------------------------
%
%"minfreq" is the minimum frequency in the ST result(Default=0)
%"maxfreq" is the maximum frequency in the ST result (Default=Nyquist)
%"samplingrate" is the time interval between samples (Default=1)
%"freqsamplingrate" is the frequency-sampling interval you desire in the ST result (Default=1)
%Passing a negative number will give the default ex. [s,t,f] = st(data,-1,-1,2,2)
%-------Outputs Returned------------------------------------------------
%
% st -a complex matrix containing the Stockwell transform.
% The rows of STOutput are the frequencies and the
% columns are the time values ie each column is
% the "local spectrum" for that point in time
% t - a vector containing the sampled times
% f - a vector containing the sampled frequencies
%--------Additional details-----------------------
% % There are several parameters immediately below that
% the user may change. They are:
%[verbose] if true prints out informational messages throughout the function.
%[removeedge] if true, removes a least squares fit parabola
% and puts a 5% hanning taper on the edges of the time series.
% This is usually a good idea.
%[analytic_signal] if the timeseries is real-valued
% this takes the analytic signal and STs it.
% This is almost always a good idea.
%[factor] the width factor of the localizing gaussian
% ie, a sinusoid of period 10 seconds has a
% gaussian window of width factor*10 seconds.
% I usually use factor=1, but sometimes factor = 3
% to get better frequency resolution.
% Copyright (c) by Bob Stockwell
% $Revision: 1.2 $ $Date: 1997/07/08 $
% This is the S transform wrapper that holds default values for the function.
TRUE = 1;
FALSE = 0;
%%% DEFAULT PARAMETERS [change these for your particular application]
verbose = TRUE;
removeedge= FALSE;
analytic_signal = FALSE;
factor = 1;
%%% END of DEFAULT PARAMETERS
%%%START OF INPUT VARIABLE CHECK
% First: make sure it is a valid time_series
% If not, return the help message
if verbose disp(' '),end % i like a line left blank
if nargin == 0
if verbose disp('No parameters inputted.'),end
st_help
t=0;,st=-1;,f=0;
return
end
% Change to column vector
if size(timeseries,2) > size(timeseries,1)
timeseries=timeseries';
end
% Make sure it is a 1-dimensional array
if size(timeseries,2) > 1
error('Please enter a *vector* of data, not matrix')
return
elseif (size(timeseries)==[1 1]) == 1
error('Please enter a *vector* of data, not a scalar')
return
end
% use defaults for input variables
if nargin == 1
minfreq = 0;
maxfreq = fix(length(timeseries)/2);
samplingrate=1;
freqsamplingrate=1;
elseif nargin==2
maxfreq = fix(length(timeseries)/2);
samplingrate=1;
freqsamplingrate=1;
[ minfreq,maxfreq,samplingrate,freqsamplingrate] = check_input(minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,timeseries);
elseif nargin==3
samplingrate=1;
freqsamplingrate=1;
[ minfreq,maxfreq,samplingrate,freqsamplingrate] = check_input(minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,timeseries);
elseif nargin==4
freqsamplingrate=1;
[ minfreq,maxfreq,samplingrate,freqsamplingrate] = check_input(minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,timeseries);
elseif nargin == 5
[ minfreq,maxfreq,samplingrate,freqsamplingrate] = check_input(minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,timeseries);
else
if verbose disp('Error in input arguments: using defaults'),end
minfreq = 0;
maxfreq = fix(length(timeseries)/2);
samplingrate=1;
freqsamplingrate=1;
end
if verbose
disp(sprintf('Minfreq = %d',minfreq))
disp(sprintf('Maxfreq = %d',maxfreq))
disp(sprintf('Sampling Rate (time domain) = %d',samplingrate))
disp(sprintf('Sampling Rate (freq. domain) = %d',freqsamplingrate))
disp(sprintf('The length of the timeseries is %d points',length(timeseries)))
disp(' ')
end
%END OF INPUT VARIABLE CHECK
% If you want to "hardwire" minfreq & maxfreq & samplingrate & freqsamplingrate do it here
% calculate the sampled time and frequency values from the two sampling rates
t = (0:length(timeseries)-1)*samplingrate;
spe_nelements =ceil((maxfreq - minfreq+1)/freqsamplingrate) ;
f = (minfreq + [0:spe_nelements-1]*freqsamplingrate)/(samplingrate*length(timeseries));
if verbose disp(sprintf('The number of frequency voices is %d',spe_nelements)),end
% The actual S Transform function is here:
st = strans(timeseries,minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,removeedge,analytic_signal,factor);
% this function is below, thus nicely encapsulated
%WRITE switch statement on nargout
% if 0 then plot amplitude spectrum
if nargout==0
if verbose disp('Plotting pseudocolor image'),end
pcolor(t,f,abs(st))
end
return
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
function st = strans(timeseries,minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,removeedge,analytic_signal,factor);
% Returns the Stockwell Transform, STOutput, of the time-series
% Code by R.G. Stockwell.
% Reference is "Localization of the Complex Spectrum: The S Transform"
% from IEEE Transactions on Signal Processing, vol. 44., number 4,
% April 1996, pages 998-1001.
%
%-------Inputs Returned------------------------------------------------
% - are all taken care of in the wrapper function above
%
%-------Outputs Returned------------------------------------------------
%
% ST -a complex matrix containing the Stockwell transform.
% The rows of STOutput are the frequencies and the
% columns are the time values
%
%
%-----------------------------------------------------------------------
% Compute the length of the data.
n=length(timeseries);
original = timeseries;
if removeedge
if verbose disp('Removing trend with polynomial fit'),end
ind = [0:n-1]';
r = polyfit(ind,timeseries,2);
fit = polyval(r,ind) ;
timeseries = timeseries - fit;
if verbose disp('Removing edges with 5% hanning taper'),end
sh_len = floor(length(timeseries)/10);
wn = hanning(sh_len);
if(sh_len==0)
sh_len=length(timeseries);
wn = 1&[1:sh_len];
end
% make sure wn is a column vector, because timeseries is
if size(wn,2) > size(wn,1)
wn=wn';
end
timeseries(1:floor(sh_len/2),1) = timeseries(1:floor(sh_len/2),1).*wn(1:floor(sh_len/2),1);
timeseries(length(timeseries)-floor(sh_len/2):n,1) = timeseries(length(timeseries)-floor(sh_len/2):n,1).*wn(sh_len-floor(sh_len/2):sh_len,1);
end
% If vector is real, do the analytic signal
if analytic_signal
if verbose disp('Calculating analytic signal (using Hilbert transform)'),end
% this version of the hilbert transform is different than hilbert.m
% This is correct!
ts_spe = fft(real(timeseries));
h = [1; 2*ones(fix((n-1)/2),1); ones(1-rem(n,2),1); zeros(fix((n-1)/2),1)];
ts_spe(:) = ts_spe.*h(:);
timeseries = ifft(ts_spe);
end
% Compute FFT's
tic;vector_fft=fft(timeseries);tim_est=toc;
vector_fft=[vector_fft,vector_fft];
tim_est = tim_est*ceil((maxfreq - minfreq+1)/freqsamplingrate) ;
if verbose disp(sprintf('Estimated time is %f',tim_est)),end
% Preallocate the STOutput matrix
st=zeros(ceil((maxfreq - minfreq+1)/freqsamplingrate),n);
% Compute the mean
% Compute S-transform value for 1 ... ceil(n/2+1)-1 frequency points
if verbose disp('Calculating S transform...'),end
if minfreq == 0
st(1,:) = mean(timeseries)*(1&[1:1:n]);
else
st(1,:)=ifft(vector_fft(minfreq+1:minfreq+n).*g_window(n,minfreq,factor));
end
%the actual calculation of the ST
% Start loop to increment the frequency point
for banana=freqsamplingrate:freqsamplingrate:(maxfreq-minfreq)
st(banana/freqsamplingrate+1,:)=ifft(vector_fft(minfreq+banana+1:minfreq+banana+n).*g_window(n,minfreq+banana,factor));
end % a fruit loop! aaaaa ha ha ha ha ha ha ha ha ha ha
% End loop to increment the frequency point
if verbose disp('Finished Calculation'),end
%%% end strans function
%------------------------------------------------------------------------
function gauss=g_window(length,freq,factor)
% Function to compute the Gaussion window for
% function Stransform. g_window is used by function
% Stransform. Programmed by Eric Tittley
%
%-----Inputs Needed--------------------------
%
% length-the length of the Gaussian window
%
% freq-the frequency at which to evaluate
% the window.
% factor- the window-width factor
%
%-----Outputs Returned--------------------------
%
% gauss-The Gaussian window
%
vector(1,:)=[0:length-1];
vector(2,:)=[-length:-1];
vector=vector.^2;
vector=vector*(-factor*2*pi^2/freq^2);
% Compute the Gaussion window
gauss=sum(exp(vector));
%-----------------------------------------------------------------------
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^%
function [ minfreq,maxfreq,samplingrate,freqsamplingrate] = check_input(minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,timeseries)
% this checks numbers, and replaces them with defaults if invalid
% if the parameters are passed as an array, put them into the appropriate variables
s = size(minfreq);
l = max(s);
if l > 1
if verbose disp('Array of inputs accepted.'),end
temp=minfreq;
minfreq = temp(1);;
if l > 1 maxfreq = temp(2);,end;
if l > 2 samplingrate = temp(3);,end;
if l > 3 freqsamplingrate = temp(4);,end;
if l > 4
if verbose disp('Ignoring extra input parameters.'),end
end;
end
if minfreq < 0 | minfreq > fix(length(timeseries)/2);
minfreq = 0;
if verbose disp('Minfreq < 0 or > Nyquist. Setting minfreq = 0.'),end
end
if maxfreq > length(timeseries)/2 | maxfreq < 0
maxfreq = fix(length(timeseries)/2);
if verbose disp(sprintf('Maxfreq < 0 or > Nyquist. Setting maxfreq = %d',maxfreq)),end
end
if minfreq > maxfreq
temporary = minfreq;
minfreq = maxfreq;
maxfreq = temporary;
clear temporary;
if verbose disp('Swapping maxfreq <=> minfreq.'),end
end
if samplingrate <0
samplingrate = abs(samplingrate);
if verbose disp('Samplingrate <0. Setting samplingrate to its absolute value.'),end
end
if freqsamplingrate < 0 % check 'what if freqsamplingrate > maxfreq - minfreq' case
freqsamplingrate = abs(freqsamplingrate);
if verbose disp('Frequency Samplingrate negative, taking absolute value'),end
end
% bloody odd how you don't end a function
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^%
function st_help
disp(' ')
disp('st() HELP COMMAND')
disp('st() returns - 1 or an error message if it fails')
disp('USAGE:: [localspectra,timevector,freqvector] = st(timeseries)')
disp('NOTE:: The function st() sets default parameters then calls the function strans()')
disp(' ')
disp('You can call strans() directly and pass the following parameters')
disp(' **** Warning! These inputs are not checked if strans() is called directly!! ****')
disp('USAGE:: localspectra = strans(timeseries,minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,removeedge,analytic_signal,factor) ')
disp(' ')
disp('Default parameters (available in st.m)')
disp('VERBOSE - prints out informational messages throughout the function.')
disp('REMOVEEDGE - removes the edge with a 5% taper, and takes')
disp('FACTOR - the width factor of the localizing gaussian')
disp(' ie, a sinusoid of period 10 seconds has a ')
disp(' gaussian window of width factor*10 seconds.')
disp(' I usually use factor=1, but sometimes factor = 3')
disp(' to get better frequency resolution.')
disp(' ')
disp('Default input variables')
disp('MINFREQ - the lowest frequency in the ST result(Default=0)')
disp('MAXFREQ - the highest frequency in the ST result (Default=nyquist')
disp('SAMPLINGRATE - the time interval between successive data points (Default = 1)')
disp('FREQSAMPLINGRATE - the number of frequencies between samples in the ST results')
% end of st_help procedure
SST.m
function [Ts] = SST(x,hlength)
% Computes the SST (Ts) of the signal x.
% INPUT
% x : Signal needed to be column vector.
% hlength: The hlength of window function.
% OUTPUT
% Ts : The SST
[xrow,xcol] = size(x);
if (xcol~=1),
error('X must be column vector');
end;
if (nargin < 1),
error('At least 1 parameter is required');
end;
if (nargin < 2),
hlength=round(xrow/5);
end;
Siglength=xrow;
hlength=hlength+1-rem(hlength,2);
ht = linspace(-0.5,0.5,hlength);ht=ht';
% Gaussian window
h = exp(-pi/0.32^2*ht.^2);
% derivative of window
dh = -2*pi/0.32^2*ht .* h; % g'
[hrow,hcol]=size(h); Lh=(hrow-1)/2;
N=xrow;
t=1:xrow;
[trow,tcol] = size(t);
tfr1= zeros (N,tcol) ;
tfr2= zeros (N,tcol) ;
tfr= zeros (round(N/2),tcol) ;
Ts= zeros (round(N/2),tcol) ;
for icol=1:tcol,
ti= t(icol); tau=-min([round(N/2)-1,Lh,ti-1]):min([round(N/2)-1,Lh,xrow-ti]);
indices= rem(N+tau,N)+1;
rSig = x(ti+tau,1);
%rSig = hilbert(real(rSig));
tfr1(indices,icol)=rSig.*conj(h(Lh+1+tau));
tfr2(indices,icol)=rSig.*conj(dh(Lh+1+tau));
end;
tfr1=fft(tfr1);
tfr2=fft(tfr2);
tfr1=tfr1(1:round(N/2),:);
tfr2=tfr2(1:round(N/2),:);
ft = 1:round(N/2);
bt = 1:N;
%%operator omega
nb = length(bt);
neta = length(ft);
va=N/hlength;
omega = zeros (round(N/2),tcol);
for b=1:nb
omega(:,b) = (ft-1)'+real(va*1i*tfr2(ft,b)/2/pi./tfr1(ft,b));
end
omega=round(omega);
for b=1:nb%time
% Reassignment step
for eta=1:neta%frequency
if abs(tfr1(eta,b))>0.0001%you can set much lower value than this.
k = omega(eta,b);
if k>=1 && k<=neta
Ts(k,b) = Ts(k,b) + tfr1(eta,b);
end
end
end
end
Ts=Ts/(xrow/2);
endSET.m
function [tfr] = SET(x,hlength)
% Synchroextracted transform
% x : Signal.
% hlength : Window length.
% tfr : Time-Frequency Representation.
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
%
% Written by Gang Yu in Shandong University at 2015.12.1.
[xrow,xcol] = size(x);
N=xrow;
if (xcol~=1),
error('X must be column vector');
end;
if (nargin < 2),
hlength=round(xrow/8);
end;
t=1:N;
ft = 1:round(N/2);
[trow,tcol] = size(t);
hlength=hlength+1-rem(hlength,2);
ht = linspace(-0.5,0.5,hlength);ht=ht';
% Gaussian window
h = exp(-pi/0.32^2*ht.^2);
% derivative of window
dh = -2*pi/0.32^2*ht .* h; % g'
[hrow,hcol]=size(h); Lh=(hrow-1)/2;
tfr1= zeros (N,tcol);
tfr2= zeros (N,tcol);
for icol=1:tcol,
ti= t(icol); tau=-min([round(N/2)-1,Lh,ti-1]):min([round(N/2)-1,Lh,xrow-ti]);
indices= rem(N+tau,N)+1;
rSig = x(ti+tau,1);
tfr1(indices,icol)=rSig.*conj(h(Lh+1+tau));
tfr2(indices,icol)=rSig.*conj(dh(Lh+1+tau));
end;
tfr1=fft(tfr1);
tfr2=fft(tfr2);
tfr1=tfr1(1:round(N/2),:);
tfr2=tfr2(1:round(N/2),:);
va=N/hlength;
omega = zeros(round(N/2),tcol);
for b=1:N
omega(:,b) = (ft-1)'+real(va*1i*tfr2(ft,b)/2/pi./tfr1(ft,b));
end
tfr=zeros(round(N/2),tcol);
for i=1:round(N/2)%frequency
for j=1:N%time
if abs(omega(i,j)-i)<0.5%default frequency resolution is 1Hz.
tfr(i,j)=tfr1(i,j);
end
end
end
tfr=tfr/(sum(h)/2);%the amplitude of tfr result has been pre-rectified.
RWT.m
function rwt = RWT(x,nvoice,wavelet,oct,scale)
% RWT -- Real Wavelet Transform
% Usage
% rwt = RWT(x,nvoice,wavelet)
% Inputs
% x signal, dyadic length n=2^J, real-valued
% nvoice number of voices/octave
% wavelet string 'Gauss', 'DerGauss','Sombrero', 'Morlet'
% Outputs
% rwt matrix n by nscale where
% nscale = nvoice .* noctave
%
% Description
% see sections 4.3.1 and 4.3.3 of Mallat's book
%
if nargin<4,
oct = 2;
scale = 4;
end
% preparation
x = ShapeAsRow(x);
n = length(x);
xhat = fft(x);
xi = [ (0: (n/2)) (((-n/2)+1):-1) ] .* (2*pi/n);
% root
omega0 = 5;
% noctave = floor(log2(n))-2;
% noctave = floor(log2(n))-1;
noctave = floor(log2(n))-oct;
nscale = nvoice .* noctave;
rwt = zeros(n,nscale);
kscale = 1;
% scale = 4;
% scale = 16;
for jo = 1:noctave,
for jv = 1:nvoice,
qscale = scale .* (2^(jv/nvoice));
omega = n .* xi ./ qscale ;
if strcmp(wavelet,'Gauss'),
window = exp(-omega.^2 ./2);
elseif strcmp(wavelet,'DerGauss'),
window = i.*omega.*exp(-omega.^2 ./2);
elseif strcmp(wavelet,'Sombrero'),
window = (omega.^2) .* exp(-omega.^2 ./2);
elseif strcmp(wavelet,'Morlet'),
window = exp(-(omega - omega0).^2 ./2) - exp(-(omega.^2 + omega0.^2)/2);
end
%( Renormalisation comme dans le bouquin )
window = window ./ sqrt(qscale);
what = window .* xhat;
w = ifft(what);
rwt(1:n,kscale) = real(w)';
kscale = kscale+1;
end
scale = scale .*2;
end
end
% kscale = 1 -> low frequencies
% the matrix rwt is ordered from low to high frequencies
% Written by Maureen Clerc and Jerome Kalifa, 1997
% clerc@cmapx.polytechnique.fr, kalifa@cmapx.polytechnique.fr
function row = ShapeAsRow(sig)
% ShapeAsRow -- Make signal a row vector
% Usage
% row = ShapeAsRow(sig)
% Inputs
% sig a row or column vector
% Outputs
% row a row vector
%
% See Also
% ShapeLike
%
row = sig(:)';
end
%
% Part of Wavelab Version 850
% Built Tue Jan 3 13:20:39 EST 2006
% This is Copyrighted Material
% For Copying permissions see COPYING.m
% Comments? e-mail wavelab@stat.stanford.edu
RelativePositionMatrix.m
function RPM=RelativePositionMatrix(x,k)
% 计算相对位置矩阵(Relative Position Matrix, RPM)
% 输入：x，一维时间序列
% k，分段聚合近似(PAA)的缩减因子
% 输出：RPM，相对位置矩阵
mu = mean(x);
delta = sqrt(var(x));
z = (x-mu)/delta;
% PAA
N = length(x);
m = ceil(N/k);
if ceil(N/k)-floor(N/k) == 0
for i = 1:m
X(i) = 1/k * sum(z(k*(i-1)+1:k*i));
end
else
for i = 1:m-1
X(i) = 1/k * sum(z(k*(i-1)+1:k*i));
end
X(m) = 1/(N-k*(m-1)) * sum(z(k*(m-1)+1:N));
end
% 计算两个时间戳之间的相对位置
M = repmat(X,m,1) - repmat(X',1,m);
%相对位移矩阵RPM
RPM = (M - min(M(:))) / (max(M(:))) - min(M(:)) * 255;
% 可视化
if nargout == 0
imagesc(RPM);
set(gcf,'Position',[300 200 450 300]); % 自行修改合适大小
xlim([0,size(RPM,1)]);
ylim([0,size(RPM,1)]);
axis square
colormap(jet); %自行修改colormap
end
endRecurrencePlot.m
function rp_mat = RecurrencePlot(x,m,tau,eta,vision)
% RecurrencePlot
% 输入数据 x
% m: 嵌入维数
% tau：时延
% eta：阈值
% vision=1 可视化
N = length(x);
T=N-(m-1)*tau;
Y=zeros(T,m);
% 相空间
for j=1:T
Y(j,:)=x(j:tau:j+(m-1)*tau);
end
%
rp_mat = ones(size(Y,1),size(Y,1)); % 初始化
for i=1:size(Y,1)
for j=i+1:size(Y,1)
d(i,j) = norm(Y(i,:)-Y(j,:));
if d(i,j) > eta
rp_mat(i,j)=0;
rp_mat(j,i)=0;
end
end
end
% 绘图
if vision==1
imshow(rp_mat,[]);%黑白
end
endMSST_new.m
function [Ts,tfr,omega2] = MSST_new(x,hlength,num)
% Computes the MSST (Ts) of the signal x.
% Expression (31)-based Algorithm.
% INPUT
% x : Signal needed to be column vector.
% hlength: The hlength of window function.
% num : iteration number.
% OUTPUT
% Ts : The SST
% tfr : The STFT
[xrow,xcol] = size(x);
if (xcol~=1),
error('X must be column vector');
end;
% if (nargin < 3),
% error('At least 3 parameter is required');
% end
%
if (nargin < 2)
hlength=round(xrow/5);
num=1;
else if (nargin < 3)
num=1;
end
end
hlength=hlength+1-rem(hlength,2);
ht = linspace(-0.5,0.5,hlength);ht=ht';
% Gaussian window
h = exp(-pi/0.32^2*ht.^2);
[hrow,~]=size(h); Lh=(hrow-1)/2;
N=xrow;
t=1:xrow;
[~,tcol] = size(t);
tfr= zeros (round(N/2),tcol) ;
omega = zeros (round(N/2),tcol-1);
omega2 = zeros (round(N/2),tcol);
Ts = zeros (round(N/2),tcol);
for icol=1:tcol,
ti= t(icol); tau=-min([round(N/2)-1,Lh,ti-1]):min([round(N/2)-1,Lh,xrow-ti]);
indices= rem(N+tau,N)+1;
rSig = x(ti+tau,1);
tfr(indices,icol)=rSig.*conj(h(Lh+1+tau));
end;
tfr=fft(tfr);
tfr=tfr(1:round(N/2),:);
for i=1:round(N/2)
omega(i,:)=diff(unwrap(angle(tfr(i,:))))*(N)/2/pi;
end
omega(:,end+1)=omega(:,end);
omega=round(omega);
[neta,nb]=size(tfr);
if num>1
for kk=1:num-1
for b=1:nb
for eta=1:neta
k = omega(eta,b);
if k>=1 && k<=neta
omega2(eta,b)=omega(k,b);
end
end
end
omega=omega2;
end
else
omega2=omega;
end
for b=1:nb%time
% Reassignment step
for eta=1:neta%frequency
if abs(tfr(eta,b))>0.0001%you can set much lower value than this.
k = omega2(eta,b);
if k>=1 && k<=neta
Ts(k,b) = Ts(k,b) + tfr(eta,b);
end
end
end
end
%tfr=tfr/(sum(h)/2);
tfr=tfr/(xrow/2);
Ts=Ts/(xrow/2);
end
%
% function [Ts_f]=SST(tfr_f,omega_f);
% [tfrm,tfrn]=size(tfr_f);
% Ts_f= zeros (tfrm,tfrn) ;
% %mx=max(max(tfr_f));
% for b=1:tfrn%time
% % Reassignment step
% for eta=1:tfrm%frequency
% %if abs(tfr_f(eta,b))>0.001*mx%you can set much lower value than this.
% k = omega_f(eta,b);
% if k>=1 && k<=tfrm
% Ts_f(k,b) = Ts_f(k,b) + tfr_f(eta,b);
% end
% %end
% end
% end
% endLMSST.m
function [Ts, IF, tfr] = LMSST(x,hlength,le)
% Local maximum synchrosqueezing transform
% Computes the SST (Ts) of the signal x.
% INPUT
% x : Signal needed to be column vector.
% hlength: The hlength of window function.
% OUTPUT
% Ts : The SST
[xrow,xcol] = size(x);
if (xcol~=1),
error('X must be column vector');
end;
if (nargin < 1),
error('At least 1 parameter is required');
end;
if (nargin < 2),
hlength=round(xrow/5);
le=round(hlength/2);
end;
if (nargin < 3),
le=round(hlength/2);
end;
%Siglength=xrow;
hlength=hlength+1-rem(hlength,2);
ht = linspace(-0.5,0.5,hlength);ht=ht';
% Gaussian window
h = exp(-pi/0.32^2*ht.^2);
% derivative of window
dh = -2*pi/0.32^2*ht .* h; % g'
[hrow,~]=size(h); Lh=(hrow-1)/2;
N=xrow;
t=1:xrow;
[~,tcol] = size(t);
tfr1= zeros (N,tcol) ;
%omega= zeros (N,tcol) ;
%tfr= zeros (round(N/2),tcol) ;
Ts= zeros (round(N/2),tcol) ;
%tfr0= ones (round(N/2),tcol) ;
IF= zeros (round(N/2),tcol) ;
for icol=1:tcol,
ti= t(icol); tau=-min([round(N/2)-1,Lh,ti-1]):min([round(N/2)-1,Lh,xrow-ti]);
indices= rem(N+tau,N)+1;
rSig = x(ti+tau,1);
%rSig = hilbert(real(rSig));
tfr1(indices,icol)=rSig.*conj(h(Lh+1+tau));
%tfr2(indices,icol)=rSig.*conj(dh(Lh+1+tau));
end;
tfr1=fft(tfr1);
%tfr2=fft(tfr2);
tfr=tfr1(1:round(N/2),:);
%tfr2=tfr2(1:round(N/2),:);
ft = 1:round(N/2);
bt = 1:N;
%%operator omega
nb = length(bt);
neta = length(ft);
%va=N/hlength;
omega = zeros (round(N/2),tcol);
%for b=1:nb
%omega(:,b) = (ft-1)'+real(va*1i*tfr2(ft,b)/2/pi./tfr1(ft,b));
%end
%omega=round(omega);
for b=1:nb%time
% Reassignment step
for eta=1:neta%frequency
%if abs(tfr1(eta,b))>0.0001%you can set much lower value than this.
[~,index]=max(abs(tfr1(max(1,eta-le):min(round(N/2),eta+le),b)));
omega(eta,b)=max(1,eta-le)+index-1-1;
%end
end
end
for b=1:nb%time
% Reassignment step
for eta=1:neta%frequency
%if abs(tfr1(eta,b))>0.0001%you can set much lower value than this.
k = omega(eta,b);
if k>=1 && k<=neta
Ts(k,b) = Ts(k,b) + tfr1(eta,b);
%SSO(k,b) = SSO(k,b) + tfr0(eta,b);
if abs(omega(eta,b)-eta)<0.5%default frequency resolution is 1Hz.
IF(eta,b)=1;
end
end
%end
end
end
Ts=Ts/(xrow/2);
endGAF.m
function GM = GAF(x,sc,vision)
% Gramian Angular Field
% 输入数据 x
% 符号控制sc选择 GASF 和 GADF
% vision=1,表示可视化，其他不可视化
norm_x = ((x-max(x))+(x-min(x)))/(max(x)-min(x)); % 归一化
% 计算
theta = acos(norm_x);
GM=zeros(length(theta));
for i = 1:length(theta)
for j = 1:length(theta)
if sc=='-'
GM(i,j) = sin(theta(i)-theta(j)); %角度减
elseif sc=='+'
GM(i,j) = cos(theta(i)+theta(j)); %角度加
end
end
end
% 绘图
if vision==1
imagesc(GM);
end
endfset.m
function [set,f,t] = fset(x,varargin)
%FSET Fourier synchroextracted transform
% SET = FSET(X) returns the Fourier synchroextracted transform of X. Each
% column of SET contains the synchroextracted spectrum of a windowed
% segment of X. The number of columns of SET is equal to the number of
% samples of the input X, which is padded with zeros on each side. Each
% spectrum is one-sided for real X and two-sided and centered for complex
% X. By default, a Kaiser window of length 256 and a beta of 10 is used.
% If X has fewer than 256 samples, then the Kaiser window has the same
% length as X.
%
% [SET,W,S] = FSET(X) returns a vector of normalized frequencies, W,
% corresponding to the rows of SET and a vector of sample numbers, S,
% corresponding to the columns of SET. Each element of S is the sample
% number of the midpoint of a windowed segment of X.
%
% [SET,F,T] = FSET(X,Fs) specifies the sampling frequency, Fs, in hertz,
% as a positive scalar. The output contains sample times, T, expressed in
% seconds and frequencies, F, expressed in hertz.
%
% [SET,F,T] = FSET(X,Ts) specifies the sampling interval, Ts, as a
% <a href="matlab:help duration">duration</a>. Ts is the time between samples of X. T has the same units
% as the duration. The units of F are in cycles/unit time of the
% duration.
%
% [...] = FSET(X,...,WINDOW) when WINDOW is a vector, divides X into
% segments of the same length as WINDOW and then multiplies each segment
% by WINDOW. If WINDOW is an integer, a Kaiser window with length WINDOW
% and a beta of 10 is used. If WINDOW is not specified, the default
% Kaiser window is used.
%
% FSET(...) with no output arguments plots the synchroextracted transform
% of the input vector x on the current figure.
%
% FSET(...,FREQLOCATION) controls where MATLAB displays the frequency
% axis on the plot. This string can be either 'xaxis' or 'yaxis'.
% Setting FREQLOCATION to 'yaxis' displays frequency on the y-axis and
% time on the x-axis. The default is 'xaxis', which displays frequency
% on the x-axis. FREQLOCATION is ignored when output arguments are
% specified.
%
% % Example 1:
% % Compute and plot the Fourier synchroextracted transform (SET) of a
% % signal that consists of two chirps.
% fs = 3000;
% t = 0:1/fs:1-1/fs;
% x1 = 2*chirp(t,500,t(end),1000);
% x2 = chirp(t,400,t(end),800);
% fset(x1+x2,fs,'yaxis')
% title('Magnitude of Fourier Synchroextracted Transform of Two Chirps')
%
% % Compute and plot the short-time Fourier transform.
% [stft,f,t] = spectrogram(x1+x2,kaiser(256,10),255,256,fs);
% figure
% h = imagesc(t,f,abs(stft));
% xlabel('Time (s)')
% ylabel('Frequency (Hz)')
% title('Magnitude of Short-time Fourier Transform of Two Chirps')
% h.Parent.YDir = 'normal';
%
% % Example 2
% % Compute the Fourier synchroextracted transform of a speech signal.
% load mtlb
% fset(mtlb,Fs,hann(256),'yaxis')
%
% See also ifsst, tfridge, spectrogram, duration.
% Copyright 2015-2018 The MathWorks, Inc.
narginchk(1,4);
nargoutchk(0,3);
if nargin > 1
[varargin{:}] = convertStringsToChars(varargin{:});
end
% Parse inputs. Fs is populated as used throughout the bulk of the code
% even if Ts is specified. fNorm specifies normalized frequencies.
[Fs,Ts,win,fNorm,freqloc] = parseInputs(x,varargin{:});
validateInputs(x,Fs,Ts,win);
% Parameters based on window size - noverlap is fixed so the transform is
% invertible.
nwin = length(win);
nfft = nwin;
noverlap = nwin-1;
% Convert to column vectors
x = signal.internal.toColIfVect(x);
win = win(:);
% Compute the output time vector (one time per sample point of the input)
if fNorm
tout = (1:length(x));
else
tout = (0:length(x)-1)/Fs;
end
% cast to enforce precision rules (we already checked that the inputs are
% numeric.
% cast to enforce precision rules
if (isa(x,'single') || isa(win,'single'))
x = single(x);
win = single(win);
Fs = single(Fs);
tout = single(tout);
else
Fs = double(Fs);
end
% Pad the signal vector x
if mod(nwin,2)
xp = [zeros((nwin-1)/2,1) ; x ; zeros((nwin-1)/2,1)];
else
xp = [zeros((nwin)/2,1) ; x ; zeros((nwin-2)/2,1)];
end
nxp = length(xp);
% Place xp into columns for the STFT
xin = getSTFTColumns(xp,nxp,nwin,noverlap,Fs);
% Compute the STFT
[sstout,fout] = computeDFT(bsxfun(@times,win,xin),nfft,Fs);
stftc = computeDFT(bsxfun(@times,dtwin(win,Fs),xin),nfft,Fs);
clear xin;
% Compute the reassignment vector
fcorr = -imag(stftc./ sstout);
fcorr(~isfinite(fcorr)) = 0;
fcorr = bsxfun(@plus,fout,fcorr);
tcorr = bsxfun(@plus,tout,zeros(size(fcorr)));
clear stftc;
% Mulitply STFT by a linear phase shift to produce the modified STFT
m = floor(nwin/2);
inds = 0:nfft-1;
ez = exp(-1i*2*pi*m*inds/nfft)';
sstout = bsxfun(@times,sstout,ez);
% Synchroextracting transform
if nargin >1
if isduration(varargin{1})
if isempty(varargin{1})
error(message('signal:fsst:EmptyDuration'));
else
Ts = varargin{1};
Fs = 1/seconds(Ts);
end
elseif ischar(varargin{1})
freqloc = validatestring(varargin{1},{'xaxis','yaxis'});
Fs=length(x);
fcorr=fcorr*Fs/2/pi;
end
else
Fs=length(x);
fcorr=fcorr*Fs/2/pi;
end
f=0:(Fs/2)/(nwin/2):round(Fs/2)-(Fs/2)/(nwin/2);
fcorr=fcorr(1:round(nwin/2),:);
[neta,nb]=size(fcorr);
set=zeros(size(fcorr));
for b=1:nb
for eta=1:neta
if abs(fcorr(eta,b)-f(eta))<0.5*(Fs/nwin)
set(eta,b)=sstout(eta,b);
end
end
end
% Reduce to one-sided spectra if the input is real, otherwise return a
% two-sided (centered) spectra.
if isreal(x)
fout = psdfreqvec('npts',nfft,'Fs',Fs,'Range','half');
sstout = sstout(1:length(fout),:);
else
% Centered spectra
sstout = centerest(sstout);
fout = centerfreq(fout);
end
% Scale fout and tout if the input is a duration object
if ~isempty(Ts)
[~,units,timeScale] = getDurationandUnits(Ts);
tout = tout*timeScale;
fout = fout/timeScale;
else
units = [];
end
if nargout == 0
if ~isempty(units)
switch units
case 'sec'
plotsst(seconds(tout),fout,set,fNorm,freqloc);
case 'min'
plotsst(minutes(tout),fout,set,fNorm,freqloc);
case 'hr'
plotsst(hours(tout),fout,set,fNorm,freqloc);
case 'day'
plotsst(days(tout),fout,set,fNorm,freqloc);
case 'year'
plotsst(years(tout),fout,set,fNorm,freqloc);
end
else
plotsst(tout,fout,set,fNorm,freqloc);
end
else
sst = sstout;
f = fout;
t = tout(:)';
end
%--------------------------------------------------------------------------
function [Fs,Ts,win,fNorm,freqloc] = parseInputs(x,varargin)
% Set defaults
Fs = [];
Ts = [];
win = kaiser(min(256,length(x)),10); %#ok<*NASGU>
fNorm = false;
freqloc = '';
% Parse optional inputs
if nargin > 1
if isduration(varargin{1})
if isempty(varargin{1})
% Throw error is empty duration object is supplied
error(message('signal:fsst:EmptyDuration'));
else
Ts = varargin{1};
Fs = 1/seconds(Ts);
end
elseif ischar(varargin{1})
freqloc = validatestring(varargin{1},{'xaxis','yaxis'});
else
if ~isempty(varargin{1})
Fs = varargin{1};
end
end
end
if isempty(Fs) && isempty(Ts)
fNorm = true;
Fs = 2*pi;
end
if nargin > 2
if ischar(varargin{2})
freqloc = validatestring(varargin{2},{'xaxis','yaxis'});
elseif ~isempty(varargin{2})
win = varargin{2};
if isscalar(win)
validateattributes(win,{'numeric'},{'positive'},'fsst','WINDOW');
win = kaiser(double(win),10);
end
end
end
if nargin > 3
freqloc = validatestring(varargin{3},{'xaxis','yaxis'});
end
%--------------------------------------------------------------------------
function validateInputs(x,Fs,Ts,win)
validateattributes(x,{'single','double'},...
{'nonsparse','finite','nonnan','vector'},'fsst','X');
validateattributes(Fs,{'numeric'},...
{'real','positive','finite','nonnan','scalar'},'fsst','Fs');
validateattributes(win,{'single','double'},...
{'real','finite','nonnegative','nonnan','vector'},'fsst','WINDOW');
if ~isempty(Ts)
dt = getDurationandUnits(Ts);
validateattributes(dt,{'numeric'},...
{'real','positive','finite','nonnan','scalar'},'fsst','Ts');
end
% Check X has at least 2 samples
if length(x) < 2
error(message('signal:fsst:MustBeMinLengthX'));
end
% Check WINDOW has at least 2 samples
if length(win) < 2
error(message('signal:fsst:MustBeMinLengthWin'));
end
% Check window length is not more than the length of the input signal.
if length(win) > length(x)
error(message('signal:fsst:WinLength'));
end
%--------------------------------------------------------------------------
function [dt,units,timeScale] = getDurationandUnits(Ts)
% This function returns the sampling interval and a format string
% The Units string is only for plotting.
tsformat = Ts.Format;
% Use first character of format string to determine correct
% duration object method.
tsformat = tsformat(1);
% Using the same time units as engunits. Units in engunits are
% not localized.
% time_units = {'secs','mins','hrs','days','years'};
switch tsformat
case 's'
dt = seconds(Ts);
units = 'sec';
timeScale = 1;
case 'm'
dt = minutes(Ts);
units = 'min';
timeScale = 1/seconds(minutes(1));
case 'h'
dt = hours(Ts);
units = 'hr';
timeScale = 1/seconds(hours(1));
case 'd'
dt = days(Ts);
units = 'day';
timeScale = 1/seconds(days(1));
case 'y'
dt = years(Ts);
units = 'year';
timeScale = 1/seconds(years(1));
end
%--------------------------------------------------------------------------
function plotsst(t,f,set,fNorm,freqloc)
% Plot the FSET in the current figure
if fNorm
%Convert to time as expected by plotTFR
t = t ./ (2*pi);
end
if isempty(freqloc)
plotOpts.freqlocation = 'xaxis';
else
plotOpts.freqlocation = freqloc;
end
plotOpts.title = getString(message('signal:fsst:titleFSST'));
plotOpts.cblbl = getString(message('signal:fsst:ColorbarLabel'));
plotOpts.cursorclbl = [getString(message('signal:fsst:ColorbarLabel')) ': '];
plotOpts.isFsnormalized = fNorm;
signalwavelet.internal.convenienceplot.plotTFR(t,f,abs(set),plotOpts);
DWVD.m
function [K,DVW,t,f] = DWVD(x,fs)
% Discrete Wigner-Ville Distribution
% x = Input Signal
% fs = sampling rate
% K = Kernel values
% DVW = DWVT
%This program is free software; you can redistribute it and/or modify
if (rem(length(x),2) == 0)
% error("Window Size must be in odd number")
x=[x;0];
end
N = length(x)-1;
X=fft(x);
X=[X(1:N/2+1);zeros(N,1);X(N/2+2:N+1)];
x=2*ifft(X);
x1 = x;
x=[zeros(N,1);x;zeros(N,1)];
z = hilbert(x);
X=zeros(2*N+1);
for k=1:2*N+1
X(:,k)=z(k+(0:2*N));
end
K = X.*conj(flipud(X));
DVW=real(fft(K([N+1:2*N+1,1:N],:)));
t = (0:2*N)/fs;
f = (0:N)'/(N+1) * fs;
f = linspace(0,f(end),2*N);
% figure
% subplot(2,2,[3,4]);
% imagesc(t,f/1e6,DVW);
% set(gca,'YDir','normal')
% xlabel('Time (s)');
% ylabel('Freqeuncy (MHz)');
% l = caxis;
% caxis([0 l(2)]);
% subplot(2,2,1);
% plot(t,real(x1));
% xlabel('Time (s)');
% ylabel('Amp');
% subplot(2,2,2);
% plot(t,max(DVW))
% xlabel('Time (S)');
% ylabel('P|x|');
end项目链接
通过网盘分享的文件：1Dto2D.zip
链接: https://pan.baidu.com/s/1MaI9OCRNN_ksBM6AcBvKzg?pwd=bpmv 提取码: bpmv
--来自百度网盘超级会员v7的分享
扫一扫
962
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
