S变换——Matlab（S变换函数，一个）

S变换——Matlab(S变换函数，一个例子)

#include 
#include 
#include 
#include 
#include 

char *Wisfile = NULL;
char *Wistemplate = “%s/.fftwis“;
#define WISLEN 8

void set_wisfile（void）
{
	char *home;

	if （Wisfile） return;
	home = getenv（“HOME“）;
	Wisfile = （char *）malloc（strlen（home） + WISLEN + 1）;
	sprintf（Wisfile Wistemplate home）;
}

/* Convert frequencies in Hz into rows of the ST given sampling rate and
length. */

int st_freq（double f int len double srate）
{
	return floor（f * len / srate + .5）;
}

static double gauss（int n int m）;

/* Stockwell transform of the real array data. The len argument is the
number of time points and it need not be a power of two. The lo and hi
arguments specify the range of frequencies to return in Hz. If they are
both zero they default to lo = 0 and hi = len / 2. The result is
returned in the complex array result which must be preallocated with
n rows and len columns where n is hi - lo + 1. For the default values of
lo and hi n is len / 2 + 1. */

void st（int len int lo int hi double *data double *result）
{
	int i k n l2;
	double s *p;
	FILE *wisdom;
	static int planlen = 0;
	static double *g;
	static fftw_plan p1 p2;
	static fftw_complex *h *H *G;

	/* Check for frequency defaults. */

	if （lo == 0 && hi == 0） {
		hi = len / 2;
	}

	/* Keep the arrays and plans around from last time since this
	is a very common case. Reallocate them if they change. */

	if （len != planlen && planlen > 0） {
		fftw_destroy_plan（p1）;
		fftw_destroy_plan（p2）;
		fftw_free（h）;
		fftw_free（H）;
		fftw_free（G）;
		free（g）;
		planlen = 0;
	}

	if （planlen == 0） {
		planlen = len;
		h = fftw_malloc（sizeof（fftw_complex） * len）;
		H = fftw_malloc（sizeof（fftw_complex） * len）;
		G = fftw_malloc（sizeof（fftw_complex） * len）;
		g = （double *）malloc（sizeof（double） * len）;

		/* Get any accumulated wisdom. */

		set_wisfile（）;
		wisdom = fopen（Wisfile “r“）;
		if （wisdom） {
			fftw_import_wisdom_from_file（wisdom）;
			fclose（wisdom）;
		}

		/* Set up the fftw plans. */

		p1 = fftw_plan_dft_1d（len h H FFTW_FORWARD FFTW_MEASURE）;
		p2 = fftw_plan_dft_1d（len G h FFTW_BACKWARD FFTW_MEASURE）;

		/* Save the wisdom. */

		wisdom = fopen（Wisfile “w“）;
		if （wisdom） {
			fftw_export_wisdom_to_file（wisdom）;
			fclose（wisdom）;
		}
	}

	/* Convert the input to complex. Also compute the mean. */

	s = 0.;
	memset（h 0 sizeof（fftw_complex） * len）;
	for （i = 0; i < len; i++） {
		h[i][0] = data[i];
		s += data[i];
	}
	s /= len;

	/* FFT. */

	fftw_execute（p1）; /* h -> H */

	/* Hilbert transform. The upper half-circle gets multiplied by
	two and the lower half-circle gets set to zero.  The real axis
	is left alone. */

	l2 = （len + 1） / 2;
	for （i = 1; i < l2; i++） {
		H[i][0] *= 2.;
		H[i][1] *= 2.;
	}
	l2 = len / 2 + 1;
	for （i = l2; i < len; i++） {
		H[i][0] = 0.;
		H[i][1] = 0.;
	}

	/* Fill in rows of the result. */

	p = result;

	/* The row for lo == 0 contains the mean. */

	n = lo;
	if （n == 0） {
		for （i = 0

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       7292  2011-08-02 22:08  S_Transform\st.c

     文件      13710  2011-08-02 22:07  S_Transform\st.m

     文件        928  2011-08-10 17:14  S_Transform\stest.asv

     文件        928  2011-08-10 17:17  S_Transform\stest.m

     目录          0  2011-08-10 17:04  S_Transform

----------- ---------  ---------- -----  ----

                22858                    5