mf_fft.hlp

是一个经济学管理应用软件 很难找的 但是经济学学生又必须用到

{smcl}
{* 31mar2005}{...}
{cmd:help mata fft()}
{hline}
{* index fft()}{...}
{* index invfft()}{...}
{* index _fft()}{...}
{* index _invfft()}{...}
{* index convolve()}{...}
{* index deconvolve()}{...}
{* index Corr()}{...}
{* index ftpad()}{...}
{* index ftwrap()}{...}
{* index ftretime()}{...}
{* index ftunwrap()}{...}
{* index ftfreqs()}{...}
{* index ftperiodogram()}{...}
{* index Fourier transform}{...}
{* index correlation}{...}

{title:Title}

{p 4 4 2}
{bf:[M-5] fft() -- Fourier transform}


{title:Syntax}

{p 8 12 2}
{it:complex vector}
{cmd:fft(}{it:numeric vector h}{cmd:)}

{p 8 12 2}
{it:numeric vector}
{cmd:invfft(}{it:numeric vector H}{cmd:)}


{p 8 12 2}
{it:void}{bind:         }
{cmd:_fft(}{it:complex vector h}{cmd:)}

{p 8 12 2}
{it:void}{bind:         }
{cmd:_invfft(}{it:complex vector H}{cmd:)}


{p 8 12 2}
{it:numeric vector}
{cmd:convolve(}{it:numeric vector r}{cmd:,}
{it:numeric vector s}{cmd:)}

{p 8 12 2}
{it:numeric vector}
{cmd:deconvolve(}{it:numeric vector r}{cmd:,}
{it:numeric vector sm}{cmd:)}


{p 8 12 2}
{it:numeric vector}
{cmd:Corr(}{it:numeric vector g}{cmd:,}
{it:numeric vector h}{cmd:,}
{it:real scalar k}{cmd:)}


{p 8 12 2}
{it:real vector}{bind:   }
{cmd:ftperiodogram(}{it:numeric vector H}{cmd:)}


{p 8 12 2}
{it:numeric vector}
{cmd:ftpad(}{it:numeric vector h}{cmd:)}

{p 8 12 2}
{it:numeric vector}
{cmd:ftwrap(}{it:numeric vector r}{cmd:,}
{it:real scalar n}{cmd:)}

{p 8 12 2}
{it:numeric vector}
{cmd:ftretime(}{it:numeric vector r}{cmd:,}
{it:numeric vector s}{cmd:)}

{p 8 12 2}
{it:numeric vector}
{cmd:ftunwrap(}{it:numeric vector H}{cmd:)}

{p 8 12 2}
{it:real vector}{bind:   }
{cmd:ftfreqs(}{it:numeric vector H}{cmd:,}
{it:real scalar delta}{cmd:)}



{title:Description}

{p 4 4 2}
{it:H}={cmd:fft(}{it:h}{cmd:)}
and 
{it:h}={cmd:invfft(}{it:H}{cmd:)} 
calculate the Fourier transform and inverse Fourier transform.
The length of {it:h} ({it:H}) must be a power of 2.

{p 4 4 2}
{cmd:_fft(}{it:h}{cmd:)}
and
{cmd:_invfft(}{it:H}{cmd:)} do the same thing, but they perform the 
calculation in place, replacing the contents of {it:h} and {it:H}.

{p 4 4 2}
{cmd:convolve(}{it:r}{cmd:,} {it:s}{cmd:)} 
returns the convolution of the signal {it:s} with the response function
{it:r}.  
{cmd:deconvolve(}{it:r}{cmd:,} {it:sm}{cmd:)}
deconvolves the smeared signal {it:sm} with the response function {it:r}, and
is thus the inverse of {cmd:convolve()}.

{p 4 4 2}
{cmd:Corr(}{it:g}{cmd:,} {it:h}{cmd:,} {it:k}{cmd:)}
returns a 2{it:k}+1 element vector containing the correlations of {it:g} and
{it:h} for lags and leads as large as {it:k}.  

{p 4 4 2}
{cmd:ftperiodogram(}{it:H}{cmd:)} returns a real vector 
containing the one-sided periodogram of {it:H}.

{p 4 4 2}
{cmd:ftpad(}{it:h}{cmd:)} returns {it:h} padded with 0s to have a length
that is a power of 2.

{p 4 4 2}
{cmd:ftwrap(}{it:r}{cmd:,} {it:n}{cmd:)} 
converts the symmetrically stored response function
{it:r} into wraparound format of length {it:n},
{it:n}>=rows({it:r})*cols({it:r}) and rows({it:r})*cols({it:r}) odd.

{p 4 4 2}
{cmd:ftunwrap(}{it:H}{cmd:)} unwraps frequency-wraparound order such as
returned by {cmd:fft()}.  You may find this useful when graphing or listing
results, but it is otherwise unnecessary.

{p 4 4 2}
{cmd:ftfreqs(}{it:H}{cmd:,} {it:delta}{cmd:)} 
returns a vector containing the frequencies associated with the elements of
{it:H}; {it:delta} is the sampling interval and is often specified as 1.


{title:Remarks}

{p 4 4 2}
Remarks are presented under the headings

	{bf:Definitions, notation, and conventions}
	{bf:Fourier transform}
	{bf:Convolution and deconvolution}
	{bf:Correlation}
	{bf:Utility routines}
	{bf:Warnings}


{title:Definitions, notation, and conventions}

{p 4 4 2}
A signal {it:h} is a row or column vector containing real or complex elements.
The length of the signal is defined as the number of elements of the vector.
It is occasionally necessary to pad a signal to a given length.  This is
done by forming a new vector equal to the original and with zeros added to the
end.

{p 4 4 2}
The Fourier transform of a signal {it:h}, typically denoted by capital letter
{it:H} of {it:h}, is stored in frequency-wraparound order.  That is, if there
are {it:n} elements in {it:H}:

			{it:H}[1]            frequency 0
			{it:H}[2]            frequency 1 
			{it:H}[3]            frequency 2
			.
			{it:H}[{it:n}/2]          frequency {it:n}/2-1
			{it:H}[{it:n}/2+1]        frequency {it:n}/2 (-{it:n}/2, aliased)
			{it:H}[{it:n}/2+2]        frequency -({it:n}/2-1)
			.
			{it:H}[{it:n}-1]          frequency -2
			{it:H}[{it:n}]            frequency -1

{p 4 4 2}
All routines expect and use this order, but see {cmd:ftunwrap()} below.

{p 4 4 2}
A response function {it:r} is a row or column vector containing
{it:m}=2{it:k}+1 real or complex elements.  {it:m} is called the duration of
the response function.  Response functions are generally stored symmetrically,
although the response function itself need not be symmetric.  The response
vector contains

		       {it:r}[1]             response at lag -{it:k}
		       {it:r}[2]             response at lag -{it:k}+1
		       .
		       {it:r}[{it:k}]             response at lag -1
		       {it:r}[{it:k}+1]           contemporaneous response 
		       {it:r}[{it:k}+2]           response at lead 1
		       {it:r}[{it:k}+3]           response at lead 2
		       .
		       {it:r}[2{it:k}+1]          response at lead {it:k}

{p 4 4 2}
Note that response functions always have odd lengths.  Response vectors are
never padded.

{p 4 4 2}
You may occasionally find it convenient to store a response vector in
"wraparound" order (similar to frequency-wraparound order), although none of
the routines here require this.  In wraparound order:

		    wrap[1]             contemporaneous response
		    wrap[2]             response at lead 1
		    wrap[3]             response at lead 2
		    .
		    wrap[{it:k}+1]           response at lead {it:k}
		    wrap[{it:k}+2]           response at lag -{it:k}
		    wrap[{it:k}+3]           response at lag -{it:k}+1
		    .
		    wrap[2{it:k}+1]          response at lag -1

{p 4 4 2}
Response vectors stored in wraparound order may be internally padded (as
opposed to merely padded) to a given length by the insertion of zeros between
wrap[{it:k}+1] and wrap[{it:k}+2].


{title:Fourier transform}

{p 4 4 2}
{cmd:fft(}{it:h}{cmd:)} returns the discrete Fourier transform of {it:h}.
{it:h} may be either real or complex, but its length must be a power of 2, 
so one typically codes {cmd:fft(ftpad(}{it:h}{cmd:))}; 
see {cmd:ftpad()}, below.  The returned result is p-conformable with {it:h}.
The calculation is performed by {cmd:_fft()}.

{p 4 4 2}
{cmd:invfft(}{it:H}{cmd:)} returns the discrete inverse Fourier transform of
{it:H}.  {it:H} may be either real or complex, but its length must be a power
of 2.  The returned result is p-conformable with {it:H}.  The calculation is
performed by {cmd:_invfft()}.

{p 4 4 2}
Note that {cmd:invfft(}{it:H}{cmd:)} may return a real or complex.  This
should be viewed as a feature, but if you wish to ensure the complex
interpretation, code {cmd:C(invfft(}{it:H}{cmd:))}.

{p 4 4 2}
{cmd:_fft(}{it:h}{cmd:)} is the built-in procedure that performs the 
fast Fourier transform in place.  {it:h} must be complex, and its length must
be a power of 2.

{p 4 4 2}
{cmd:_invfft(}{it:H}{cmd:)} is the built-in procedure that performs the
inverse fast Fourier transform in place.  {it:H} must be complex, and its
length must be a power of 2.


{title:Convolution and deconvolution}

{p 4 4 2}
{cmd:convolve(}{it:r}{cmd:,} {it:s}{cmd:)} 
returns the convolution of the signal {it:s} with the response function
{it:r}.  Calculation is performed by taking the {cmd:fft()} of the elements,
multiplying, and using {cmd:invfft()} to transform the results back.
Nevertheless, it is not necessary that the length of {it:s} be a power of 2.
{cmd:convolve()} handles all paddings necessary, including paddings to {it:s}
required to prevent the result from being contaminated by erroneous wrapping
around of {it:s}.  Although one thinks of the convolution operator as
being commutative, {cmd:convolve()} is not commutative since required
zero-padding of the response and signal differ.

{p 4 4 2}
If {it:n} is the length of the signal and 2{it:k}+1 is the length of the
response function, the returned result has length {it:n}+2{it:k}.  The first 
{it:k} elements are the convoluted signal before the true signal begins, and
the last {it:k} elements are the convoluted signal after the true signal ends.
See {cmd:ftretime()}, below.  In any case, you may only be interested in the
elements {cmd:convolve()[|}{it:k}+1{cmd:\}{it:n}-{it:k}{cmd:|]}, 
the part contemporaneous with {it:s}.

{p 4 4 2}
The returned vector is a row vector if {it:s} is a row vector and a column
vector otherwise.  The result is guaranteed to be real if both {it:r} and
{it:s} are real; the result may be complex or real, otherwise.

{p 4 4 2}
It is not required that the response function be shorter than the signal,
although this will typically be the case.


{p 4 4 2}
{cmd:deconvolve(}{it:r}{cmd:,} {it:sm}{cmd:)}
deconvolves the smeared signal {it:sm} with the response function {it:r}, and
is thus the inverse of {cmd:convolve()}.  In particular,

		{cmd:deconvolve(}{it:r}{cmd:, convolve(}{it:r}{cmd:,}{it:s}{cmd:))} == {it:s}        (up to roundoff error)