#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?		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?		H[i][0]?*=?2.;
		H[i][1]?*=?2.;
	}
	l2?=?len?/?2?+?1;
	for?（i?=?l2;?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

?屬性????????????大小?????日期????時(shí)間???名稱
-----------?---------??----------?-----??----

?????文件???????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