S变换介绍(附代码)
1、S变换
作为小波变换和短时傅里叶变换的继承和发展, S 变换采用高斯窗函数且窗宽与频率的倒数成正比,免去了窗函数的选择和改善了窗宽固定的缺陷,并且时频表示中各频率分量的相位谱与原始信号保持直接的联系,使其在 PQD 分析中可以采用更多的特征量,同时, S 变换提取的特征量对噪声不敏感。因此,在电能电能质量扰动、轴承故障诊断领域运用广泛。适用于非平稳信号的时频分析方法,其定义为:
2、S变换matlab代码
unction [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
代码引用来源IEEE Transactions on Signal Processing, vol. 44., number 4, April 1996, pages 998-1001.
S变换介绍(附代码)相关推荐
- php 运营商授权,PHP判断手机号运营商(详细介绍附代码)
道理很简单,知道手机号规则 进行正则判断就可以 移动:134.135.136.137.138.139.150.151.157(TD).158.159.187.188 联通:130.131.132.15 ...
- 【自动驾驶】30.c++实现基于eigen实现欧拉角(RPY), 旋转矩阵, 旋转向量, 四元数之间的变换(附代码)
矩阵的使用可参考系列博客:点击此处 原文链接:基于eigen实现欧拉角(RPY), 旋转矩阵, 旋转向量, 四元数之间的变换. 也可以参考另一篇博客:eigen 中四元数.欧拉角.旋转矩阵.旋转向量. ...
- 一个模型搞定多个CTR业务!阿里STAR网络介绍(附代码实现)
今天为大家带来阿里巴巴2021年的一篇文章:<One Model to Serve All: Star Topology Adaptive Recommender for Multi-Domai ...
- STM32学习笔记:读写内部Flash(介绍+附代码)
一.介绍 首先我们需要了解一个内存映射: stm32的flash地址起始于0x0800 0000,结束地址是0x0800 0000加上芯片实际的flash大小,不同的芯片flash大小不同. RAM起 ...
- 数字图像处理Matlab-图像的滤波处理与图像空间变换(附代码)
目录 1.Objectives: 2.Experiment Content: 3.Experiment Principle: 4.Experiment Steps Result and Conlusi ...
- STM32读写内部Flash(介绍+附代码)
概述 内部flash读写详解 一.介绍 首先我们需要了解一个内存映射: stm32的flash地址起始于0x0800 0000,结束地址是0x0800 0000加上芯片实际的flash大小,不同的芯片 ...
- 详细介绍附代码:使用jquery,和php文件构建一个简单的在线聊天室,通过ip显示googlemap
最近学习了关于使用最为流行的jquery发送请求,在实践中以最为简单的聊天室作为测验的辅助工具,对相关网页开发有一个初步的认识,希望大家能够一起学习进步. 首先介绍一下相关文件信息和功能 ...
- php 判断号码运营商,PHP根据手机号判断运营商(详细介绍附代码)
道理很简单,知道手机号规则 进行正则判断就可以 移动:134.135.136.137.138.139.150.151.157(TD).158.159.187.188 联通:130.131.132.15 ...
- php判断运营商,PHP根据手机号判断运营商(详细介绍附代码)
道理很简单,知道手机号规则 进行正则判断就可以 移动:134.135.136.137.138.139.150.151.157(TD).158.159.187.188 联通:130.131.132.15 ...
- C# List排序简介及四种方法介绍-附代码
有时用户需要按某项排序,但是查询结果以List格式存储,我们当然可以自己编写一个快速排序的方法进行排序,但是还有多个选择,并且可能比你写的短.效率也不差,那不如在恰当的时候选择其他方法对List进行排 ...
最新文章
- 使用 Fries 创建性感的 Android 风格移动应用界面
- 统计与分布之高斯分布
- 华三交换机升级的ipe文件_弱电工程工业以太网交换机电源故障总结
- vim 如何将特定范围行注释掉,以及在末尾添加注释
- 重温强化学习之无模型学习方法:时间差分方法
- java父进程获取子进程异常_如何在perl的父进程中获取死亡的子进程的PID?
- PHP自学3——在html的table标签中显示用户提交表单
- AIOps智能化数据体系的构建及在字节跳动的实践
- mac mysql常用命令
- java h5 调用摄像头_基于百度AI使用H5实现调用摄像头进行人脸注册、人脸搜索功能(Java)...
- sprint演示会议
- 科技论文之Introduction部分写作
- 计算机网络位置设置工作组,工作组设置【处置步骤】
- iOS 地图坐标系转换
- 互联网公司的完整开发流程是怎样的?
- hdu 5296 Annoying problem (LCA)
- 基于JavaWeb医疗管理系统的开发与实现
- Android国外学习资源汇总
- 少年,请多一些开疆拓土的勇气——写给在C和C++间犹豫的学生
- 【pandas】教程:1-处理什么样的数据