S(Stockwell)变换的Matlab代码

S变换简介

S变换，又称为Stockwell变换，由R. G. Stockwell于1996年提出。具体的定义如下：

S变换在傅里叶域的表示形式为：

离散的S变换为：

S变换克服了短时傅里叶变换固定窗函数宽度的缺陷，采用了一个随频率变化的高斯窗函数。它的窗函数宽度与频率的倒数成正比，高频时用窄窗，低频用宽窗。所以具有多分辨率分析的思想，也可以看作相位校正的小波变换。

Stockwell版S变换程序

function [st, t, f] = st(timeseries, minfreq, maxfreq, samplingrate, freqsamplingrate)
% Returns the Stockwell Transform of the timeseries.
% Code by Robert Glenn 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 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 = FALSE;
removeedge = FALSE;
analytic_signal = TRUE;
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));
%%
lamda=15;%5  15
p=1;%0.5
vector(1,:)=[0:length-1];
vector(2,:)=[-length:-1];
vector=vector.^2;
vector=vector*(-factor*lamda^2*pi^2/(freq.^(2*p)));
gauss=sum(lamda*freq.^p*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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389

Dash版S变换程序

function ST=stran(h)

% Compute S-Transform without for loops

%%% Coded by Kalyan S. Dash %%%
%%% IIT Bhubaneswar, India %%%

[~,N]=size(h); % h is a 1xN one-dimensional series

nhaf=fix(N/2);

odvn=1;

if nhaf*2==N
    odvn=0;
end

f=[0:nhaf -nhaf+1-odvn:-1]/N;

Hft=fft(h);

%Compute all frequency domain Gaussians as one matrix

invfk=[1./f(2:nhaf+1)]';

W=2*pi*repmat(f,nhaf,1).*repmat(invfk,1,N);

G=exp((-W.^2)/2); %Gaussian in freq domain

% End of frequency domain Gaussian computation

% Compute Toeplitz matrix with the shifted fft(h)

HW=toeplitz(Hft(1:nhaf+1)',Hft);

% Exclude the first row, corresponding to zero frequency

HW=[HW(2:nhaf+1,:)];

% Compute Stockwell Transform

ST=ifft(HW.*G,[],2); %Compute voice

%Add the zero freq row

st0=mean(h)*ones(1,N);

ST=[st0;ST];

end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50

博主自己编写的S变换

function wcoefs = myst(t,Sig,freqlow,freqhigh,alpha)
%==================%
% Stockwell transform ver 1.0 (2021-09-24)
% Copyright (c) by Jinshun Shen,Xidian University
%jsshen@stu.xidian.edu.cn
%This program can not be used for commercialization without the authorization of its author
%======输入======%
%t:时间序列
% Sig: 输入信号
%freqlow,freqhigh:频率范围
%alpha:频率分辨率
%======输出======%
% wcoefs: ST变换计算结果

if (nargin <= 2)

error('At least 2 parameters required');

end

if (nargin < 5)
    alpha=0.05;
elseif (nargin < 4)
    freqhigh=log(length(t));
elseif (nargin < 3)
    freqlow=0;
end

if size(t,1)>1%t是列向量,则转置
    t=t';
end
if size(Sig,1)>1%Sig是列向量,则转置
    Sig=Sig';
end
% 信号的长度
SigLen = length(Sig);
% 时间的长度
TimeLen = length(t);
dt=t(2)-t(1);
fre=freqlow:alpha:freqhigh;%产生频率范围
% time=1:TimeLen;
if SigLen > TimeLen%信号长度大于时间长度，时间补0
    t=[t  zeros(1,SigLen-TimeLen)];
end
if SigLen < TimeLen%信号长度小于于时间长度，信号补0
    Sig=[Sig  zeros(1,TimeLen-SigLen)];
end
% 总频率数量
nLevel=length(fre);
% 分配计算结果的存储单元
wcoefs = zeros(nLevel,TimeLen); 
temp=zeros(1,TimeLen);
for m = 1:nLevel
% 计算各频率上的ST系数
% 提取频率参数
f= fre(m);
    for n=1:TimeLen
        %计算高斯窗函数
        Gauss_st=(f/(sqrt(2*pi)))*exp(-0.5*(n*dt-t).^2*(f)^2).*exp(-1i*2*pi*f*t);
        temp(n)=trapz(t,Sig.*Gauss_st);%效率低，精度稍高
    end
wcoefs(m,:)=temp;
end
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64

R. G. Stockwell, L. Mansinha, and R. Lowe, “Localization of the
complex spectrum: the S transform,” IEEE transactions on signal
processing, vol. 44, no. 4, pp. 998–1001, 1996.