声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 3206|回复: 4

[小波] 同步压缩小波变换(sst)

[复制链接]
发表于 2016-12-4 16:58 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
近段时间一直在看同步压缩小波变换,但是效果甚微,本人matlab基础很差,现在附上sst的程序,希望大神指点一下,怎么确定程序中的几个量(s,gamma,sigma,ft,bt)谢谢了。
function [Vs Ts Rs] = sst(s,gamma,sigma,ft,bt)
% computes the STFT Vs, the SST Ts and the reassigned STFT Rs of the signal s.
% Uses a Gaussian window with SD sigma. Computes the TF transforms at
% frequencies ft and discrete times bt.
%
% Input:
%   s: signal
%   gamma: threshold
%   sigma: window parameter
%   ft: index of frequencies where to compute the transforms.
%       default: 1:N/2 (N=length(s))
%   bt: time index where to compute the transforms
%       default: 1:N
% checking length of signal
n = length(s);
nv = log2(n);
if mod(nv,1)~=0
    warning('The signal is not a power of two, truncation to the next power');
    s = s(1:2^floor(nv));
end
n = length(s);
s = s(:);
% Optional parameters
if nargin<5
   ft = 1:n/2;
   bt = 1:n;
end
% Padding
sleft = flipud(s(2:n/2+1));
sright = flipud(s(end-n/2:end-1));
x = [sleft; s; sright];
clear xleft xright;
%% STFT and operators tau and omega
nb = length(bt);
neta = length(ft);
t = linspace(-0.5,0.5,n);t=t';
g = exp(-pi/sigma^2*t.^2);
gp = -2*pi/sigma^2*t .* g; % g'
Vs = zeros(neta,nb);
Ts = zeros(neta,nb);
Rs = zeros(neta,nb);
%% Direct reassignment
df = ft(2)-ft(1);
db = bt(2)-bt(1);
for b=1:nb
    % Computing Vs, omega and tau in t=b
    tmp = fft(x(bt(b):bt(b)+n-1).*g)/n/sigma;
    vtmp = tmp(ft);
    tmp = fft(x(bt(b):bt(b)+n-1).*gp)/n/sigma;
    omegtmp = (ft-1)'-real(tmp(ft)/2/1i/pi./vtmp);
   
    tmp = fft(x(bt(b):bt(b)+n-1).*g)/1/sigma;
    tmp = [tmp(2)-tmp(1); (tmp(3:end)-tmp(1:end-2))/2; tmp(end)-tmp(end-1)];
    tautmp = bt(b)-1+real(tmp(ft)/2/pi/1i ./vtmp);
   
   
    % Storing values
    Vs(:,b) = vtmp;
   
    % Reassignment step
for eta=1:neta
        if abs(vtmp(eta))>gamma
            k = 1+round((omegtmp(eta)-ft(1)+1)/df);
            if k>=1 && k<=neta
                Ts(k,b) = Ts(k,b) + vtmp(eta)*exp(1i*pi*(ft(eta)-1));
                l = 1+round((tautmp(eta)-bt(1)+1)/db);
                if l>=1 && l<=nb
                    Rs(k,l) = Rs(k,l) + abs(vtmp(eta));
                end
            end
        end
    end
end
Ts = df*sigma*Ts;
Rs = Rs * df * db;
end

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

发表于 2016-12-7 08:35 | 显示全部楼层
你应该把运行中出现的问题发一下  不然大家都要再去运行
 楼主| 发表于 2016-12-11 12:45 | 显示全部楼层
Triste 发表于 2016-12-7 08:35
你应该把运行中出现的问题发一下  不然大家都要再去运行

请问大神,你有sst中的test.m文件吗,关于信号重构的,真的谢谢了。

点评

抱歉 这个真没有  详情 回复 发表于 2016-12-12 08:31
 楼主| 发表于 2016-12-11 12:47 | 显示全部楼层
2581428947 发表于 2016-12-11 12:45
请问大神,你有sst中的test.m文件吗,关于信号重构的,真的谢谢了。

比如信号为y=100sin(100pi*t)+30sin(250pi*t)+20sin(350pi*t)
发表于 2016-12-12 08:31 | 显示全部楼层
2581428947 发表于 2016-12-11 12:45
请问大神,你有sst中的test.m文件吗,关于信号重构的,真的谢谢了。

抱歉  这个真没有
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-25 00:11 , Processed in 0.090054 second(s), 20 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表