S变换在特征提取中的使用

S变换

S变换采用高斯窗函数且窗宽与频率的倒数成正比，免去了窗函数的选择和改善了窗宽固定的缺陷，并且时频表示中各频率分量的相位谱与原始信号保持直接的联系，S变换具有良好的时频特性，适合用S变换对信号的一些时频与特征进行提取。

可以看出S变换不同于短时傅里叶变换之处在于高斯窗的高度和宽度随频率变化，这样就克服了短时傅里叶变化窗口高度和宽度固定的缺陷。S变化的离散变现形式：

S变换的结果

S变换的结果是一个二维矩阵，列对应采样点时间，行对应频率值，矩阵元素是一个复数，可以从这个元素里面获取幅值和相位信息。需要说明的是此时的频率值需要根据采样频率和采样数进行转化才能和实际的频率对应，转换方式与离散傅里叶转化方式相同，每个点的实际频率是(fs/N)*n，fs为采样频率，N为采样数，n为第n个点的序号。S变化的结果可以用三维立体图、二维等高线以及灰度图等直观表示，或者截取某一剖面用二维曲线表示。

S变化提取特征

国内现阶段在电能质量监测和评价方面根据S变化判断电压突降、突升、谐波以及震荡等电能扰动类型。通过s变化对不同电压扰动进行时频域特征提取很流行，合理选择分类器进行训练分类，最后对电压扰动进行检测与判断。下面就以一组电力信号为例，进行S变换特征提取以及结果简要分析。

给定一个电力信号序列，对该序列进行S变换分析，获取一个二维矩阵。首先用三维图画出各个分量，如下所示：

由于电流序列的基频是60Hz,如果对该信号进行谐波分析，就需要对谐波含量进行提取，从上图可以看出信号中含有较大的直流成分。由于谐波中大多是奇次谐波，偶次谐波通常只会发生在电路谐振的条件下。下面就从二维矩阵中抽取一行采样点与幅值信息，给出频率为60Hz(基波)时的幅值曲线。

再根据二维矩阵获取该序列信号的频谱图：

通过以上图可以得知该信号中的谐波分量的大致分布。接下来可以根据实际情况选择S变化的形态域、频域、时频域特征作为样本特征，选择合适的分类器进行训练。

下面给出一段开源的S变换仿真代码：


function [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 messageif verbose disp(' '),end  % i like a line left blankif nargin == 0 if verbose disp('No parameters inputted.'),endst_helpt=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) > 1error('Please enter a *vector* of data, not matrix')return
elseif (size(timeseries)==[1 1]) == 1error('Please enter a *vector* of data, not a scalar')return
end% use defaults for input variablesif nargin == 1minfreq = 0;maxfreq = fix(length(timeseries)/2);samplingrate=1;freqsamplingrate=1;
elseif nargin==2maxfreq = 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'),endminfreq = 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'),endpcolor(t,f,abs(st))
endreturn
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 removeedgeif verbose disp('Removing trend with polynomial fit'),endind = [0:n-1]';r = polyfit(ind,timeseries,2);fit = polyval(r,ind) ;timeseries = timeseries - fit;if verbose disp('Removing edges with 5% hanning taper'),endsh_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 isif size(wn,2) > size(wn,1)wn=wn';   endtimeseries(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_signalif 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 == 0st(1,:) = mean(timeseries)*(1&[1:1:n]);
elsest(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.'),endtemp=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.'),endend;end      if minfreq < 0 | minfreq > fix(length(timeseries)/2);minfreq = 0;if verbose disp('Minfreq < 0 or > Nyquist. Setting minfreq = 0.'),endendif maxfreq > length(timeseries)/2  | maxfreq < 0 maxfreq = fix(length(timeseries)/2);if verbose disp(sprintf('Maxfreq < 0 or > Nyquist. Setting maxfreq = %d',maxfreq)),endendif minfreq > maxfreq temporary = minfreq;minfreq = maxfreq;maxfreq = temporary;clear temporary;if verbose disp('Swapping maxfreq <=> minfreq.'),endendif samplingrate <0samplingrate = abs(samplingrate);if verbose disp('Samplingrate <0. Setting samplingrate to its absolute value.'),endendif freqsamplingrate < 0   % check 'what if freqsamplingrate > maxfreq - minfreq' casefreqsamplingrate = abs(freqsamplingrate);if verbose disp('Frequency Samplingrate negative, taking absolute value'),endend% bloody odd how you don't end a function%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^%
function st_helpdisp(' ')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

S变换在特征提取中的使用相关推荐

计算机视觉技术在图像特征提取中的应用研究,基于图像特征提取的图像融合研究...
基于图像特征提取的图像融合研究 [摘要]:视觉信息是人类从自然界中获取信息的最主要手段,图像信息是一种主观性很强的重要信息表达形式,也是最难由计算机认知.处理与实现的信息之一.而图像特征提取作为计算机 ...
Tech Tuesday 22-5-31 小波分析在信号特征提取中的应用（1）
Tech Tuesday 22-5-31 文献阅读报告(1)--小波分析在信号特征提取中的应用原文下载链接点击此处 About Tech Tuesday LeBron James有Taco Tues ...
基于python的图像Gabor变换及特征提取
基于python的图像Gabor变换及特征提取 1.前言 2. "Gabor帮主"简介 3."Gabor帮主"大招之图像变换 3."Gabor帮主&q ...
《矩阵与变换》教学中的几个“务必”
<矩阵与变换>教学中的几个"务必" 摘要:新课程改革,<矩阵与变换>进入高中课堂.无论是对于富有教学经验的老教师,还是即将走向高中讲台的新教师,<矩阵 ...
hough变换检测直线 matlab,Matlab实现Hough变换检测图像中的直线
Hough变换的原理: 将图像从图像空间变换至参数空间,变换公式如下: 变换以后,图像空间与参数空间存在以下关系: 图像空间中的一点在参数空间是一条曲线,而图像空间共线的各点对应于参数空间交于一点的各 ...
小波包变换/能量特征提取/结果图绘制-python代码
1. 小波外部包下载要下载两个包: PyWavelets和Matplotlib(要运行PyWavelets的所有测试,您还需要安装 Matplotlib软件包.) 下载方法: pip install ...
计算机视觉技术在图像特征提取中的应用研究,彩色图像特征提取研究(一)
彩色图像特征提取研究(一) 彩色图像特征提取研究徐红霞摘要:以普通的彩色图像为例,介绍了对彩色图像特征提取的原理.其具体过程分为原图像的预处理.图像信息分析.图像的特征提取,然后用MATLAB实现 ...
FFT变换频谱图中幅值的设置方法
按照上篇博文所画出来的频谱图中,原信号的每个频率是准确地找出来了,但是各个频率点所对应的的幅值可不是原信号中真正的幅值,因为在进行DFT(FFT)变换的时候,已经把幅值改变了,要想让频谱图的纵坐标显示 ...
FFT变换频谱图中频率刻度的设置方法
看到matlab中关于fft变换的几行代码,总想把它们几行语句搞清楚,看了许多,还是有些搞不清楚,可能需要更多的知识才能把它们彻底搞懂吧. 先来看一个简单的画频谱图的代码吧: clear all fs ...
dct图像压缩c语言实现,DCT变换在图像压缩中的实现
小白拙见,希望理解不对的地方大家多多指教! 对于各种信号,都可以说它是由多个振幅与频率不同的正弦或者余弦函数组成的.并且一个信号通常由一个直流信号DC(幅值保持不变的信号)和多个交流信号AC(幅值以某 ...

S变换在特征提取中的使用

S变换

可以看出S变换不同于短时傅里叶变换之处在于高斯窗的高度和宽度随频率变化，这样就克服了短时傅里叶变化窗口高度和宽度固定的缺陷。S变化的离散变现形式：

S变换的结果

S变化提取特征

给定一个电力信号序列，对该序列进行S变换分析，获取一个二维矩阵。首先用三维图画出各个分量，如下所示：

再根据二维矩阵获取该序列信号的频谱图：

通过以上图可以得知该信号中的谐波分量的大致分布。接下来可以根据实际情况选择S变化的形态域、频域、时频域特征作为样本特征，选择合适的分类器进行训练。

下面给出一段开源的S变换仿真代码：

S变换在特征提取中的使用相关推荐

最新文章

热门文章