一种新的时频分析方法—同步提取变换SET

hehehala

“笔者提出了一种新的时频分析方法，取名“同步提取变换”(Synchroextracting Transform，SET)，也是为实现ITFA做的一种努力。单从名字上看，这个方法与同步压缩变换，同出一门，但是二者原理南辕东辙”
同步压缩也好（sst）,其他相关的推广,例如chui还有其他人的二阶同步压缩或者基于凸优化的同步压缩重点和难点就是在瞬时频率的提取，而你是在人家的基础上做了个小小的改变：对提取后的瞬时频率的赋值方式改变了一小下，操作上更简便，聚集性更好点，这能说南辕北辙？至于你的博士论文里说你的变换能避免不确定性原理，建议好好重读daubechies等人的工作，看看在提取瞬时频率的时候能不能规避不确定性原理。
劝你一句，搞懂了再吹。

impulse

hehehala 发表于 2017-8-25 08:40
“笔者提出了一种新的时频分析方法，取名“同步提取变换”(Synchroextracting Transform，SET)，也是为实现 ...

积极讨论使大家都进步！

g-123

hehehala 发表于 2017-8-25 08:40

南辕东辙南辕东辙

Jack315

精华，全是精华……得好好学习。

yaozirel

楼主太牛了，我要下载下来学习学习了，顺便问一下二阶系统的阶跃响应用该方法处理可行否？

yangqijiang

楼主说的对，加油

Jing666

好想法，好东西，顶起来

kevin19821

楼主，想问下，和HHT比较起来，优点在哪里？另外你在国内的期刊上发表过文章没？

补充内容 (2018-1-9 14:14):
楼主，程序中的tau=-min([round(N/2)-1,Lh,ti-1]):min([round(N/2)-1,Lh,xrow-ti])没看懂，能解释下吗

kevin19821

clc
clear all
close all
fs=1200;
dt=1/fs;
n=2000;
t=0:dt:(n-1)*dt;
s1=2*sin(2*pi*50*t);
s2=sin(2*pi*150*t);
s=s1+s2;
s(1000:1500)=0;
fre=(fs/2)/(n/2):(fs/2)/(n/2):(fs/2);
figure(1)
plot(t,s);
title('原始信号','FontSize',12);
ylabel('幅值(Hz)','FontSize',12);
xlabel('时间(t/s)','FontSize',12);
ss=s';
[srow,scol]=size(ss);
hlength=srow/5;
hlength=hlength+1-rem(hlength,2);%%%%保证hlength为奇数，后面的程序要用到
ht=linspace(-0.5,0.5,hlength);ht=ht';
h=exp(-pi/0.32^2*ht.^2);   %%%高斯窗
dh=-2*pi/0.32^2*ht .* h;
figure(2)
subplot(211)
plot(ht,h);title('GAUSSIN WINDOW');
subplot(212)
plot(ht,dh);title('微分高斯窗');
[wrow wcol]=size(h);
Lh=(wrow-1)/2;
[trow,tcol] = size(t);
tfr1=zeros (srow,srow) ;
tfr2=zeros (srow,srow) ;
tfr=zeros (round(srow/2),tcol) ;
Ts=zeros (round(srow/2),tcol) ;
t=1:srow;

N=srow;
s=s';
t=1:N;
for icol=1:tcol,
ti= t(icol); tau=-min([round(N/2)-1,Lh,ti-1]):min([round(N/2)-1,Lh,srow-ti]);
indices= rem(N+tau,N)+1;
rSig = s(ti+tau,1);
tfr1(indices,icol)=rSig.*conj(h(Lh+1+tau));
tfr2(indices,icol)=rSig.*conj(dh(Lh+1+tau));
end   %%%此循环相当于固定窗，进行信号的移动，实现STFT.
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;
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));   %%%此处没有看懂，论文中的，怎么用(ft-1)'代替了吗？并且ft=1:N,怎么用这样的一个数和后面的相加。另这里的omega应该是瞬时频率，取实部（real)主要是做什么用的？
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/(srow/2);
% end
f=(fs/2)/(n/2):(fs/2)/(n/2):(fs/2);
figure(3)
plot(f,abs(Ts));

kevin19821

C:\Users\Administrator\Desktop

kevin19821

楼主，能否解释下下面两个疑问：
（1）omega(:,b) = (ft-1)'+real(va*1i*tfr2(ft,b)/2/pi./tfr1(ft,b));   %%%此处没有看懂，论文中的，怎么用(ft-1)'代替了吗？并且ft=1:N,怎么用这样的一个数和后面的相加。另这里的omega应该是瞬时频率，取实部（real)主要是做什么用的？
（2）   Ts(k,b) = Ts(k,b) + tfr1(eta,b);%%%同步压缩变换的公式为。，这句和这个公式看起来不对应呀。

g-123

麻烦版主帮忙删除Synchroextracting Transform.pdf 这个附件，涉及到了杂志侵权问题，不好意思。

迷你猫

占个位置留个言，虽然现在还看不太明白，万一火了呢！

g-123

kevin19821 发表于 2018-1-12 09:35
楼主，能否解释下下面两个疑问：
（1）omega(:,b) = (ft-1)'+real(va*1i*tfr2(ft,b)/2/pi./tfr1(ft,b));   ...

ft-1就是文中的w频率，之所以-1是因为fft会有个1的漂移。
Omega本身是复数，只有取实部才能使用。
积分本身就是加法，这个公式就是积分。

zswseu

这个比较牛，不知道算法有没有坚实的数学理论支撑。