哪位帮看下同步压缩变换的程序

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);

hi凌寒

楼主现在在做什么，这个代码弄懂了吗

guojuntao0192

omega(:,b) = (ft-1)'+real(va*1i*tfr2(ft,b)/2/pi./tfr1(ft,b));   %这里为什么要乘以va?