scilab/modules/signal_processing/help/en_US/transforms/fft.xml

   1 <?xml version="1.0" encoding="UTF-8"?>
   2 <!--
   3  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   4  * Copyright (C) 1997   - INRIA
   5  * Copyright (C) 2012 - Serge Steer - INRIA
   6  *
   7  * This file must be used under the terms of the CeCILL.
   8  * This source file is licensed as described in the file COPYING, which
   9  * you should have received as part of this distribution.  The terms
  10  * are also available at
  11  * http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
  12  *
  13  -->
  14 <refentry xmlns="http://docbook.org/ns/docbook"
  15           xmlns:xlink="http://www.w3.org/1999/xlink"
  16           xmlns:svg="http://www.w3.org/2000/svg"
  17           xmlns:mml="http://www.w3.org/1998/Math/MathML"
  18           xmlns:db="http://docbook.org/ns/docbook"
  19           xmlns:scilab="http://www.scilab.org" xml:lang="en" xml:id="fft">
  20     <refnamediv>
  21         <refname>fft</refname>
  22         <refpurpose>fast Fourier transform.</refpurpose>
  23     </refnamediv>
  24     <refnamediv xml:id="ifft">
  25         <refname>ifft</refname>
  26         <refpurpose>fast Fourier transform.</refpurpose>
  27     </refnamediv>
  28     <refsynopsisdiv>
  29         <title>Calling Sequence</title>
  30         <synopsis>X=fft(A [,sign] [,option])
  31             X=fft(A,sign,selection  [,option])
  32             X=fft(A,sign,dims,incr [,option] )
  33         </synopsis>
  34     </refsynopsisdiv>
  35     <refsection>
  36         <title>Arguments</title>
  37         <variablelist>
  38             <varlistentry>
  39                 <term>A</term>
  40                 <listitem>
  41                     <para>a real or complex vector or real or complex array
  42                         (vector, matrix or N-D array).
  43                     </para>
  44                 </listitem>
  45             </varlistentry>
  46             <varlistentry>
  47                 <term>X</term>
  48                 <listitem>
  49                     a real or complex array with same shape as <literal>A</literal>.
  50                 </listitem>
  51             </varlistentry>
  52             <varlistentry>
  53                 <term>sign</term>
  54                 <listitem>
  55                     an integer. with possible values <literal>1</literal> or
  56                     <literal>-1</literal>. Select direct or inverse
  57                     transform. The default value is <literal>-1</literal>
  58                     (direct transform).
  59                 </listitem>
  60             </varlistentry>
  61             <varlistentry>
  62                 <term>option</term>
  63                 <listitem>
  64                     a character string. with possible values
  65                     <literal>"symmetric"</literal> or
  66                     <literal>"nonsymmetric"</literal>. Indicates if
  67                     <literal>A</literal> is symmetric or not. If this argument
  68                     is omitted the algorithm automatically determines if
  69                     <literal>A</literal> is symmetric or not. See the
  70                     Description part for details.
  71                 </listitem>
  72             </varlistentry>
  73             <varlistentry>
  74                 <term>selection</term>
  75                 <listitem>
  76                     a vector containing index on <literal>A</literal> array
  77                     dimensions.  See the Description part for details.
  78                 </listitem>
  79             </varlistentry>
  80
  81             <varlistentry>
  82                 <term>dims</term>
  83                 <listitem>
  84                     a vector of positive numbers with integer values, or a
  85                     vector of positive integers.  See the Description part for details.
  86                     <para>
  87                         Each element must be a divisor
  88                         of the total number of elements of <literal>A</literal>.
  89                     </para>
  90                     <para>
  91                         The product of the elements must be less than the total
  92                         number of elements of <literal>A</literal>.
  93                     </para>
  94                 </listitem>
  95             </varlistentry>
  96             <varlistentry>
  97                 <term>incr</term>
  98                 <listitem>
  99                     a vector of positive numbers with integer values, or a
 100                     vector of positive integers.  See the Description part for
 101                     details.
 102                     <para>
 103                         <literal>incr</literal> must have the same number of
 104                         elements than <literal>dims</literal>.
 105                     </para>
 106                     <para>
 107                         Each element must be a divisor of the total number of
 108                         elements of <literal>A</literal>.
 109                     </para>
 110                     <para>
 111                         The <literal>incr</literal> elements must be in strictly
 112                         increasing order.
 113                     </para>
 114                 </listitem>
 115             </varlistentry>
 116         </variablelist>
 117     </refsection>
 118     <refsection>
 119         <title>Description</title> This function realizes direct or
 120         inverse 1-D or N-D Discrete Fourier Transforms.
 121         <variablelist>
 122             <varlistentry>
 123                 <term>Short syntax </term>
 124                 <listitem>
 125                     <variablelist>
 126                         <varlistentry>
 127                             <term>direct</term>
 128                             <listitem>
 129                                 <literal>X=fft(A,-1 [,option])</literal> or
 130                                 <literal>X=fft(A [,option])</literal> gives a direct
 131                                 transform.
 132                                 <variablelist>
 133                                     <varlistentry>
 134                                         <term>single variate</term>
 135                                         <listitem>
 136                                             <para>
 137                                                 If <literal>A</literal> is a vector a
 138                                                 single variate direct FFT is computed that
 139                                                 is:
 140                                             </para>
 141                                             <para>
 142                                                 <latex>
 143                                                     $x(k) = \sum_{m=1}^n {a(m)*e^{-\frac{2i*\pi}{n}(m-1) (k-1)}$
 144                                                 </latex>
 145                                             </para>
 146                                             <para>
 147                                                 (the <literal>-1</literal> argument refers
 148                                                 to the sign of the exponent..., NOT to
 149                                                 "inverse"),
 150                                             </para>
 151                                         </listitem>
 152                                     </varlistentry>
 153
 154                                     <varlistentry>
 155                                         <term>multivariate</term>
 156                                         <listitem>
 157                                             <para>
 158                                                 If <literal>A</literal> is a matrix or
 159                                                 a multidimensional array a multivariate direct
 160                                                 FFT is performed.
 161                                             </para>
 162                                         </listitem>
 163                                     </varlistentry>
 164                                 </variablelist>
 165
 166                             </listitem>
 167                         </varlistentry>
 168                         <varlistentry>
 169                             <term>inverse</term>
 170                             <listitem>
 171                                 <para>
 172                                     <literal>X=fft(A,1)</literal> or
 173                                     <literal>X=ifft(A)</literal>performs the inverse
 174                                     normalized transform, such
 175                                     that<literal>A==ifft(fft(A))</literal>.
 176                                 </para>
 177                                 <variablelist>
 178                                     <varlistentry>
 179                                         <term>single variate</term>
 180                                         <listitem>
 181                                             If <literal>A</literal> is a vector a single
 182                                             variate inverse FFT is computed
 183                                             <para>
 184                                                 <latex>
 185                                                     $x(k) = \sum_{m=1}^n
 186                                                     {a(m)*e^{+\frac{2i*\pi}{n} (m-1) (k-1)}$
 187                                                 </latex>
 188                                             </para>
 189                                         </listitem>
 190                                     </varlistentry>
 191                                     <varlistentry>
 192                                         <term>multivariate</term>
 193                                         <listitem>
 194                                             <para>
 195                                                 If <literal>a</literal> is a matrix or or
 196                                                 a multidimensional array a multivariate inverse
 197                                                 FFT is performed.
 198                                             </para>
 199                                         </listitem>
 200                                     </varlistentry>
 201                                 </variablelist>
 202                             </listitem>
 203                         </varlistentry>
 204                     </variablelist>
 205                 </listitem>
 206             </varlistentry>
 207
 208             <varlistentry>
 209                 <term>Long syntax for FFT along specified dimensions</term>
 210                 <listitem>
 211                     <itemizedlist>
 212                         <listitem>
 213                             <para>
 214                                 <literal>X=fft(A,sign,selection [,option])</literal>
 215                                 allows to perform efficiently all direct or inverse
 216                                 fft of the "slices" of <literal>A</literal> along
 217                                 selected dimensions.
 218                             </para>
 219                             <para>
 220                                 For example, if <literal>A</literal> is a 3-D array
 221                                 <literal>X=fft(A,-1,2)</literal> is equivalent to:
 222                             </para>
 223                             <programlisting role=""><![CDATA[
 224 for i1=1:size(A,1)
 225   for i3=1:size(A,3)
 226     X(i1,:,i3)=fft(A(i1,:,i3),-1);
 227   end
 228 end
 229               ]]></programlisting>
 230                             <para>
 231                                 and <literal>X=fft(A,-1,[1 3])</literal> is equivalent to:
 232                             </para>
 233                             <programlisting role=""><![CDATA[
 234 for i2=1:size(A,2)
 235   X(:,i2,:)=fft(A(:,i2,:),-1);
 236 end
 237               ]]></programlisting>
 238                         </listitem>
 239                         <listitem>
 240                             <para>
 241                                 <literal>X=fft(A,sign,dims,incr [,option])</literal> is
 242                                 a previous syntax that also allows to perform all direct or
 243                                 inverse fft of the slices of <literal>A</literal> along
 244                                 selected dimensions.
 245                             </para>
 246                             <para>
 247                                 For example, if <literal>A</literal> is an array with
 248                                 <literal>n1*n2*n3</literal> elements
 249                                 <literal>X=fft(A,-1,n1,1)</literal> is equivalent to
 250                                 <literal>X=fft(matrix(A,[n1,n2,n3]),-1,1)</literal>.
 251                                 and <literal>X=fft(A,-1,[n1 n3],[1 n1*n2])</literal>
 252                                 is equivalent to
 253                                 <literal>X=fft(matrix(A,[n1,n2,n3]),-1,[1,3])</literal>.
 254                             </para>
 255                         </listitem>
 256                     </itemizedlist>
 257                 </listitem>
 258             </varlistentry>
 259             <varlistentry>
 260                 <term>Using option argument</term> This argument can be used
 261                 to inform the fft algorithm about the symmetry of
 262                 <literal>A</literal> or of all its "slices".  An N-D array
 263                 <literal>B</literal> with dimensions <literal>n1</literal>,
 264                 ..., <literal>np</literal> is conjugate symmetric for the fft
 265                 if and only if <literal>B==conj(B([1 n1:-1:2],[1
 266                     n2:-1:2],...,[1 np:-1:2]))
 267                 </literal>
 268                 .In such a case the
 269                 result <literal>X</literal> is real and an efficient specific
 270                 algorithm can be used.
 271                 <listitem>
 272                     <itemizedlist>
 273                         <listitem>
 274                             <term>"symmetric"</term> that value causes fft to treat
 275                             <literal>A</literal> or all its "slices" conjugate
 276                             symmetric. This option is useful to avoid automatic
 277                             determination of symmetry or if <literal>A</literal> or
 278                             all its "slices" are not exactly symmetric because of
 279                             round-off errors.
 280                         </listitem>
 281                         <listitem>
 282                             <term>"nonsymmetric"</term> that value causes fft not to
 283                             take care of symmetry. This option is useful to avoid
 284                             automatic determination of symmetry.
 285                         </listitem>
 286                         <listitem>
 287                             <term>unspecified</term> If the option is omitted the
 288                             fft algorithm automatically checks for exact symmetry.
 289                         </listitem>
 290
 291                     </itemizedlist>
 292                 </listitem>
 293             </varlistentry>
 294             <varlistentry>
 295                 <term>Optimizing fft</term>
 296                 <listitem>
 297                     <para>
 298                         Remark: fftw function automatically stores his last
 299                         parameters in memory to re-use it in a second time. This
 300                         improves greatly the time computation when consecutives
 301                         calls (with same parameters) are performed.
 302                     </para>
 303                     <para>
 304                         It is possible to go further in fft optimization using
 305                         <link linkend="get_fftw_wisdom">get_fftw_wisdom</link>, <link
 306              linkend="set_fftw_wisdom">set_fftw_wisdom</link> functions.
 307                     </para>
 308                 </listitem>
 309             </varlistentry>
 310         </variablelist>
 311     </refsection>
 312     <refsection>
 313         <title>Algorithms</title>
 314         <para>
 315             This function uses the  <ulink url="http://www.fftw.org/">fftw3</ulink> library.
 316         </para>
 317     </refsection>
 318     <refsection>
 319         <title>Examples</title>
 320         <para>1-D fft</para>
 321         <programlisting role="example"><![CDATA[
 322 //Frequency components of a signal
 323 //----------------------------------
 324 // build a noised signal sampled at 1000hz  containing  pure frequencies
 325 // at 50 and 70 Hz
 326 sample_rate=1000;
 327 t = 0:1/sample_rate:0.6;
 328 N=size(t,'*'); //number of samples
 329 s=sin(2*%pi*50*t)+sin(2*%pi*70*t+%pi/4)+grand(1,N,'nor',0,1);
 330
 331 y=fft(s);
 332
 333 //s is real so the fft response is conjugate symmetric and we retain only the first N/2 points
 334 f=sample_rate*(0:(N/2))/N; //associated frequency vector
 335 n=size(f,'*')
 336 clf()
 337 plot(f,abs(y(1:n)))
 338      ]]></programlisting>
 339         <para>2-D fft</para>
 340         <programlisting role="example"><![CDATA[
 341 ----------------------------------
 342 A = zeros(256,256);
 343 A(5:24,13:17) = 1;
 344 X = fftshift(fft(A));
 345 set(gcf(),"color_map",jetcolormap(128));
 346 clf;grayplot(0:255,0:255,abs(X)')
 347      ]]></programlisting>
 348         <para>mupliple fft</para>
 349         <programlisting role="example"><![CDATA[
 350 //simple case, 3 1-D fft at a time
 351 N=2048;
 352 t=linspace(0,10,2048);
 353 A=[2*sin(2*%pi*3*t)+ sin(2*%pi*3.5*t)
 354   10*sin(2*%pi*8*t)
 355    sin(2*%pi*0.5*t)+4*sin(2*%pi*0.8*t)];
 356 X=fft(A,-1,2);
 357
 358 fs=1/(t(2)-t(1));
 359 f=fs*(0:(N/2))/N; //associated frequency vector
 360 clf;plot(f(1:100)',abs(X(:,1:100))')
 361 legend(["3 and 3.5 Hz","8 Hz","0.5 and 0.8 Hz"],"in_upper_left")
 362
 363 // 45  3-D fft at a time
 364 Dims=[5 4 9 5 6];
 365 A=matrix(rand(1,prod(Dims)),Dims);
 366
 367 y=fft(A,-1,[2 4 5]);
 368
 369 //equivalent (but less efficient code)
 370 y1=zeros(A);
 371 for i1=1:Dims(1)
 372   for i3=1:Dims(3)
 373     ind=list(i1,:,i3,:,:);
 374     y1(ind(:))=fft(A(ind(:)),-1);
 375   end
 376 end
 377    ]]></programlisting>
 378
 379         <programlisting role="example"><![CDATA[
 380 //Using explicit formula for  1-D discrete Fourier transform
 381 //------------------------------------------------
 382 function xf=DFT(x,flag);
 383   n=size(x,'*');
 384   //Compute the n by n Fourier matrix
 385   if flag==1 then,//backward transformation
 386     am=exp(2*%pi*%i*(0:n-1)'*(0:n-1)/n);
 387   else //forward transformation
 388     am=exp(-2*%pi*%i*(0:n-1)'*(0:n-1)/n);
 389   end
 390   xf=am*matrix(x,n,1);//dft
 391   xf=matrix(xf,size(x));//reshape
 392   if flag==1 then,xf=xf/n;end
 393 endfunction
 394
 395 //Comparison with the fast Fourier algorithm
 396 a=rand(1,1000);
 397 norm(DFT(a,1) - fft(a,1))
 398 norm(DFT(a,-1) - fft(a,-1))
 399
 400 timer();DFT(a,-1);timer()
 401 timer();fft(a,-1);timer()
 402      ]]></programlisting>
 403     </refsection>
 404     <refsection role="see also">
 405         <title>See Also</title>
 406         <simplelist type="inline">
 407             <member>
 408                 <link linkend="corr">corr</link>
 409             </member>
 410             <member>
 411                 <link linkend="fftw_flags">fftw_flags</link>
 412             </member>
 413             <member>
 414                 <link linkend="get_fftw_wisdom">get_fftw_wisdom</link>
 415             </member>
 416             <member>
 417                 <link linkend="set_fftw_wisdom">set_fftw_wisdom</link>
 418             </member>
 419             <member>
 420                 <link linkend="fftw_forget_wisdom">fftw_forget_wisdom</link>
 421             </member>
 422         </simplelist>
 423     </refsection>
 424     <refsection>
 425         <title>Bibliography</title>
 426         <para>
 427             Matteo Frigo and Steven G. Johnson, "FFTW Documentation" <ulink url="http://www.fftw.org/#documentation">http://www.fftw.org/#documentation</ulink>
 428         </para>
 429     </refsection>
 430 </refentry>