| //+------------------------------------------------------------------+
|
| //| fasttransforms.mqh |
|
| //| Copyright 2003-2022 Sergey Bochkanov (ALGLIB project) |
|
| //| Copyright 2012-2023, MetaQuotes Ltd. |
|
| //| https://www.mql5.com |
|
| //+------------------------------------------------------------------+
|
| //| Implementation of ALGLIB library in MetaQuotes Language 5 |
|
| //| |
|
| //| The features of the library include: |
|
| //| - Linear algebra (direct algorithms, EVD, SVD) |
|
| //| - Solving systems of linear and non-linear equations |
|
| //| - Interpolation |
|
| //| - Optimization |
|
| //| - FFT (Fast Fourier Transform) |
|
| //| - Numerical integration |
|
| //| - Linear and nonlinear least-squares fitting |
|
| //| - Ordinary differential equations |
|
| //| - Computation of special functions |
|
| //| - Descriptive statistics and hypothesis testing |
|
| //| - Data analysis - classification, regression |
|
| //| - Implementing linear algebra algorithms, interpolation, etc. |
|
| //| in high-precision arithmetic (using MPFR) |
|
| //| |
|
| //| This file is free software; you can redistribute it and/or |
|
| //| modify it under the terms of the GNU General Public License as |
|
| //| published by the Free Software Foundation (www.fsf.org); either |
|
| //| version 2 of the License, or (at your option) any later version. |
|
| //| |
|
| //| This program is distributed in the hope that it will be useful, |
|
| //| but WITHOUT ANY WARRANTY; without even the implied warranty of |
|
| //| MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
|
| //| GNU General Public License for more details. |
|
| //+------------------------------------------------------------------+
|
| #include "alglibinternal.mqh"
|
| //+------------------------------------------------------------------+
|
| //| Fast Fourier Transform |
|
| //+------------------------------------------------------------------+
|
| class CFastFourierTransform
|
| {
|
| public:
|
| static void FFTC1D(complex &a[],const int n);
|
| static void FFTC1DInv(complex &a[],const int n);
|
| static void FFTR1D(double &a[],const int n,complex &f[]);
|
| static void FFTR1DInv(complex &f[],const int n,double &a[]);
|
| static void FFTR1DInternalEven(double &a[],const int n,double &buf[],CFtPlan &plan);
|
| static void FFTR1DInvInternalEven(double &a[],const int n,double &buf[],CFtPlan &plan);
|
| static void FFTC1D(CRowComplex &a,const int n);
|
| static void FFTC1DInv(CRowComplex &a,const int n);
|
| static void FFTR1D(CRowDouble &a,const int n,CRowComplex &f);
|
| static void FFTR1DInv(CRowComplex &f,const int n,CRowDouble &a);
|
| static void FFTR1DInternalEven(CRowDouble &a,const int n,CRowDouble &buf,CFtPlan &plan);
|
| static void FFTR1DInvInternalEven(CRowDouble &a,const int n,CRowDouble &buf,CFtPlan &plan);
|
| };
|
| //+------------------------------------------------------------------+
|
| //| 1-dimensional complex FFT. |
|
| //| Array size N may be arbitrary number (composite or prime). |
|
| //| Composite N's are handled with cache-oblivious variation of a |
|
| //| Cooley-Tukey algorithm. Small prime-factors are transformed using|
|
| //| hard coded codelets (similar to FFTW codelets, but without |
|
| //| low-level optimization), large prime-factors are handled with |
|
| //| Bluestein's algorithm. |
|
| //| Fastests transforms are for smooth N's (prime factors are 2, 3, |
|
| //| 5 only), most fast for powers of 2. When N have prime factors |
|
| //| larger than these, but orders of magnitude smaller than N, |
|
| //| computations will be about 4 times slower than for nearby highly |
|
| //| composite N's. When N itself is prime, speed will be 6 times |
|
| //| lower. |
|
| //| Algorithm has O(N*logN) complexity for any N (composite or |
|
| //| prime). |
|
| //| INPUT PARAMETERS |
|
| //| A - array[0..N-1] - complex function to be transformed |
|
| //| N - problem size |
|
| //| OUTPUT PARAMETERS |
|
| //| A - DFT of a input array, array[0..N-1] |
|
| //| A_out[j] = SUM(A_in[k]*exp(-2*pi*sqrt(-1)*j*k/N), |
|
| //| k = 0..N-1) |
|
| //+------------------------------------------------------------------+
|
| void CFastFourierTransform::FFTC1D(complex &a[],const int n)
|
| {
|
| //--- create a variable
|
| int i=0;
|
| //--- create array
|
| CRowDouble buf;
|
| //--- object of class
|
| CFtPlan plan;
|
| //--- check
|
| if(!CAp::Assert(n>0,__FUNCTION__+": incorrect N!"))
|
| return;
|
| //--- check
|
| if(!CAp::Assert(CAp::Len(a)>=n,__FUNCTION__+": Length(A)<N!"))
|
| return;
|
| //--- check
|
| if(!CAp::Assert(CApServ::IsFiniteComplexVector(a,n),__FUNCTION__+": A contains infinite or NAN values!"))
|
| return;
|
| //--- Special case: N=1,FFT is just identity transform.
|
| //--- After this block we assume that N is strictly greater than 1.
|
| if(n==1)
|
| return;
|
| //--- convert input array to the more convinient format
|
| buf.Resize(2*n);
|
| for(i=0; i<n; i++)
|
| {
|
| buf.Set(2*i+0,a[i].real);
|
| buf.Set(2*i+1,a[i].imag);
|
| }
|
| //--- Generate plan and execute it.
|
| //--- Plan is a combination of a successive factorizations of N and
|
| //--- precomputed data. It is much like a FFTW plan,but is not stored
|
| //--- between subroutine calls and is much simpler.
|
| CFtBase::FtComplexFFTPlan(n,1,plan);
|
| CFtBase::FtApplyPlan(plan,buf,0,1);
|
| //--- result
|
| for(i=0; i<n; i++)
|
| {
|
| a[i].real=buf[2*i+0];
|
| a[i].imag=buf[2*i+1];
|
| }
|
| }
|
| //+------------------------------------------------------------------+
|
| //| |
|
| //+------------------------------------------------------------------+
|
| void CFastFourierTransform::FFTC1D(CRowComplex &a,const int n)
|
| {
|
| //--- create a variable
|
| int i=0;
|
| //--- create array
|
| CRowDouble buf;
|
| //--- object of class
|
| CFtPlan plan;
|
| //--- check
|
| if(!CAp::Assert(n>0,__FUNCTION__+": incorrect N!"))
|
| return;
|
| //--- check
|
| if(!CAp::Assert(CAp::Len(a)>=n,__FUNCTION__+": Length(A)<N!"))
|
| return;
|
| //--- check
|
| if(!CAp::Assert(CApServ::IsFiniteComplexVector(a,n),__FUNCTION__+": A contains infinite or NAN values!"))
|
| return;
|
| //--- Special case: N=1,FFT is just identity transform.
|
| //--- After this block we assume that N is strictly greater than 1.
|
| if(n==1)
|
| return;
|
| //--- convert input array to the more convinient format
|
| buf.Resize(2*n);
|
| for(i=0; i<n; i++)
|
| {
|
| buf.Set(2*i+0,a[i].real);
|
| buf.Set(2*i+1,a[i].imag);
|
| }
|
| //--- Generate plan and execute it.
|
| //--- Plan is a combination of a successive factorizations of N and
|
| //--- precomputed data. It is much like a FFTW plan,but is not stored
|
| //--- between subroutine calls and is much simpler.
|
| CFtBase::FtComplexFFTPlan(n,1,plan);
|
| CFtBase::FtApplyPlan(plan,buf,0,1);
|
| //--- result
|
| for(i=0; i<n; i++)
|
| {
|
| a.SetRe(i,buf[2*i+0]);
|
| a.SetIm(i,buf[2*i+1]);
|
| }
|
| }
|
| //+------------------------------------------------------------------+
|
| //| 1-dimensional complex inverse FFT. |
|
| //| Array size N may be arbitrary number (composite or prime). |
|
| //| Algorithm has O(N*logN) complexity for any N (composite or prime)|
|
| //| See FFTC1D() description for more information about algorithm |
|
| //| performance. |
|
| //| INPUT PARAMETERS |
|
| //| A - array[0..N-1] - complex array to be transformed |
|
| //| N - problem size |
|
| //| OUTPUT PARAMETERS |
|
| //| A - inverse DFT of a input array, array[0..N-1] |
|
| //| A_out[j] = SUM(A_in[k]/N*exp(+2*pi*sqrt(-1)*j*k/N), |
|
| //| k = 0..N-1) |
|
| //+------------------------------------------------------------------+
|
| void CFastFourierTransform::FFTC1DInv(complex &a[],const int n)
|
| {
|
| //--- create a variable
|
| int i=0;
|
| //--- check
|
| if(!CAp::Assert(n>0,__FUNCTION__+": incorrect N!"))
|
| return;
|
| //--- check
|
| if(!CAp::Assert(CAp::Len(a)>=n,__FUNCTION__+": Length(A)<N!"))
|
| return;
|
| //--- check
|
| if(!CAp::Assert(CApServ::IsFiniteComplexVector(a,n),__FUNCTION__+": A contains infinite or NAN values!"))
|
| return;
|
| //--- Inverse DFT can be expressed in terms of the DFT as
|
| //--- invfft(x)=fft(x')'/N
|
| //--- here x' means conj(x).
|
| for(i=0; i<=n-1; i++)
|
| a[i].imag=-a[i].imag;
|
| //--- function call
|
| FFTC1D(a,n);
|
| //--- change values
|
| for(i=0; i<=n-1; i++)
|
| {
|
| a[i].real=a[i].real/n;
|
| a[i].imag=-(a[i].imag/n);
|
| }
|
| }
|
| //+------------------------------------------------------------------+
|
| //| |
|
| //+------------------------------------------------------------------+
|
| void CFastFourierTransform::FFTC1DInv(CRowComplex &a,const int n)
|
| {
|
| //--- create a variable
|
| int i=0;
|
| //--- check
|
| if(!CAp::Assert(n>0,__FUNCTION__+": incorrect N!"))
|
| return;
|
| //--- check
|
| if(!CAp::Assert(CAp::Len(a)>=n,__FUNCTION__+": Length(A)<N!"))
|
| return;
|
| //--- check
|
| if(!CAp::Assert(CApServ::IsFiniteComplexVector(a,n),__FUNCTION__+": A contains infinite or NAN values!"))
|
| return;
|
| //--- Inverse DFT can be expressed in terms of the DFT as
|
| //--- invfft(x)=fft(x')'/N
|
| //--- here x' means conj(x).
|
| for(i=0; i<n; i++)
|
| a.SetIm(i,-a[i].imag);
|
| //--- function call
|
| FFTC1D(a,n);
|
| //--- change values
|
| for(i=0; i<=n-1; i++)
|
| {
|
| a.SetRe(i,a[i].real/n);
|
| a.SetIm(i,-(a[i].imag/n));
|
| }
|
| }
|
| //+------------------------------------------------------------------+
|
| //| 1-dimensional real FFT. |
|
| //| Algorithm has O(N*logN) complexity for any N (composite or |
|
| //| prime). |
|
| //| INPUT PARAMETERS |
|
| //| A - array[0..N-1] - real function to be transformed |
|
| //| N - problem size |
|
| //| OUTPUT PARAMETERS |
|
| //| F - DFT of a input array, array[0..N-1] |
|
| //| F[j] = SUM(A[k]*exp(-2*pi*sqrt(-1)*j*k/N), |
|
| //| k = 0..N-1) |
|
| //| NOTE: |
|
| //| F[] satisfies symmetry property F[k] = conj(F[N-k]), so just |
|
| //| one half of array is usually needed. But for convinience |
|
| //| subroutine returns full complex array (with frequencies above |
|
| //| N/2), so its result may be used by other FFT-related subroutines.|
|
| //+------------------------------------------------------------------+
|
| void CFastFourierTransform::FFTR1D(double &a[],const int n,complex &f[])
|
| {
|
| //--- create variables
|
| int i=0;
|
| int n2=0;
|
| int idx=0;
|
| complex hn=0;
|
| complex hmnc=0;
|
| complex v=0;
|
| int i_=0;
|
| //--- create array
|
| CRowDouble buf;
|
| //--- object of class
|
| CFtPlan plan;
|
| //--- check
|
| if(!CAp::Assert(n>0,__FUNCTION__+": incorrect N!"))
|
| return;
|
| //--- check
|
| if(!CAp::Assert(CAp::Len(a)>=n,__FUNCTION__+": Length(A)<N!"))
|
| return;
|
| //--- check
|
| if(!CAp::Assert(CApServ::IsFiniteVector(a,n),__FUNCTION__+": A contains infinite or NAN values!"))
|
| return;
|
| //--- Special cases:
|
| //--- * N=1,FFT is just identity transform.
|
| //--- * N=2,FFT is simple too
|
| //--- After this block we assume that N is strictly greater than 2
|
| if(n==1)
|
| {
|
| //--- allocation
|
| ArrayResize(f,1);
|
| f[0]=a[0];
|
| //--- exit the function
|
| return;
|
| }
|
| //--- check
|
| if(n==2)
|
| {
|
| //--- allocation
|
| ArrayResize(f,2);
|
| f[0].real=a[0]+a[1];
|
| f[0].imag=0;
|
| f[1].real=a[0]-a[1];
|
| f[1].imag=0;
|
| //--- exit the function
|
| return;
|
| }
|
| //--- Choose between odd-size and even-size FFTs
|
| if(n%2==0)
|
| {
|
| //--- even-size real FFT,use reduction to the complex task
|
| n2=n/2;
|
| //--- allocation
|
| buf=a;
|
| buf.Resize(n);
|
| //--- function call
|
| CFtBase::FtComplexFFTPlan(n2,1,plan);
|
| CFtBase::FtApplyPlan(plan,buf,0,1);
|
| //--- allocation
|
| ArrayResize(f,n);
|
| //--- calculation
|
| for(i=0; i<=n2; i++)
|
| {
|
| idx=2*(i%n2);
|
| hn.real=buf[idx+0];
|
| hn.imag=buf[idx+1];
|
| idx=2*((n2-i)%n2);
|
| hmnc.real=buf[idx+0];
|
| hmnc.imag=-buf[idx+1];
|
| v.real=-MathSin(-(2*M_PI*i/n));
|
| v.imag=MathCos(-(2*M_PI*i/n));
|
| f[i]=hn+hmnc-v*(hn-hmnc);
|
| f[i].real=0.5*f[i].real;
|
| f[i].imag=0.5*f[i].imag;
|
| }
|
| //--- copy
|
| for(i=n2+1; i<=n-1; i++)
|
| f[i]=CMath::Conj(f[n-i]);
|
| }
|
| else
|
| {
|
| //--- use complex FFT
|
| ArrayResize(f,n);
|
| //--- copy
|
| for(i=0; i<=n-1; i++)
|
| f[i]=a[i];
|
| //--- function call
|
| FFTC1D(f,n);
|
| }
|
| }
|
| //+------------------------------------------------------------------+
|
| //| |
|
| //+------------------------------------------------------------------+
|
| void CFastFourierTransform::FFTR1D(CRowDouble &a,const int n,CRowComplex &f)
|
| {
|
| //--- create variables
|
| int i=0;
|
| int n2=0;
|
| int idx=0;
|
| complex hn=0;
|
| complex hmnc=0;
|
| complex v=0;
|
| int i_=0;
|
| //--- create array
|
| CRowDouble buf;
|
| //--- object of class
|
| CFtPlan plan;
|
| //--- check
|
| if(!CAp::Assert(n>0,__FUNCTION__+": incorrect N!"))
|
| return;
|
| //--- check
|
| if(!CAp::Assert(CAp::Len(a)>=n,__FUNCTION__+": Length(A)<N!"))
|
| return;
|
| //--- check
|
| if(!CAp::Assert(CApServ::IsFiniteVector(a,n),__FUNCTION__+": A contains infinite or NAN values!"))
|
| return;
|
| //--- Special cases:
|
| //--- * N=1,FFT is just identity transform.
|
| //--- * N=2,FFT is simple too
|
| //--- After this block we assume that N is strictly greater than 2
|
| if(n==1)
|
| {
|
| //--- allocation
|
| f.Resize(1);
|
| f.Set(0,a[0]);
|
| //--- exit the function
|
| return;
|
| }
|
| //--- check
|
| if(n==2)
|
| {
|
| //--- allocation
|
| f.Resize(2);
|
| f.Set(0,a[0]+a[1]);
|
| f.Set(1,a[0]-a[1]);
|
| //--- exit the function
|
| return;
|
| }
|
| //--- Choose between odd-size and even-size FFTs
|
| if(n%2==0)
|
| {
|
| //--- even-size real FFT,use reduction to the complex task
|
| n2=n/2;
|
| //--- allocation
|
| buf=a;
|
| buf.Resize(n);
|
| //--- function call
|
| CFtBase::FtComplexFFTPlan(n2,1,plan);
|
| CFtBase::FtApplyPlan(plan,buf,0,1);
|
| //--- allocation
|
| f.Resize(n);
|
| //--- calculation
|
| for(i=0; i<=n2; i++)
|
| {
|
| idx=2*(i%n2);
|
| hn.real=buf[idx+0];
|
| hn.imag=buf[idx+1];
|
| idx=2*((n2-i)%n2);
|
| hmnc.real=buf[idx+0];
|
| hmnc.imag=-buf[idx+1];
|
| v.real=-MathSin(-(2*M_PI*i/n));
|
| v.imag=MathCos(-(2*M_PI*i/n));
|
| f.Set(i,hn+hmnc-v*(hn-hmnc));
|
| f.SetRe(i,0.5*f[i].real);
|
| f.SetIm(i,0.5*f[i].imag);
|
| }
|
| //--- copy
|
| for(i=n2+1; i<n; i++)
|
| f.Set(i,CMath::Conj(f[n-i]));
|
| }
|
| else
|
| {
|
| //--- use complex FFT
|
| f.Resize(n);
|
| //--- copy
|
| for(i=0; i<n; i++)
|
| f.Set(i,a[i]);
|
| //--- function call
|
| FFTC1D(f,n);
|
| }
|
| }
|
| //+------------------------------------------------------------------+
|
| //| 1-dimensional real inverse FFT. |
|
| //| Algorithm has O(N*logN) complexity for any N (composite or |
|
| //| prime). |
|
| //| INPUT PARAMETERS |
|
| //| F - array[0..floor(N/2)] - frequencies from forward real |
|
| //| FFT |
|
| //| N - problem size |
|
| //| OUTPUT PARAMETERS |
|
| //| A - inverse DFT of a input array, array[0..N-1] |
|
| //| NOTE: |
|
| //| F[] should satisfy symmetry property F[k] = conj(F[N-k]), |
|
| //| so just one half of frequencies array is needed - elements from 0|
|
| //| to floor(N/2). F[0] is ALWAYS real. If N is even F[floor(N/2)] is|
|
| //| real too. If N is odd, then F[floor(N/2)] has no special |
|
| //| properties. |
|
| //| Relying on properties noted above, FFTR1DInv subroutine uses only|
|
| //| elements from 0th to floor(N/2)-th. It ignores imaginary part of |
|
| //| F[0], and in case N is even it ignores imaginary part of |
|
| //| F[floor(N/2)] too. |
|
| //| When you call this function using full arguments list - |
|
| //| "FFTR1DInv(F,N,A)" |
|
| //| - you can pass either either frequencies array with N elements or|
|
| //| reduced array with roughly N/2 elements - subroutine will |
|
| //| successfully transform both. |
|
| //| If you call this function using reduced arguments list - |
|
| //| "FFTR1DInv(F,A)" - you must pass FULL array with N elements |
|
| //| (although higher N/2 are still not used) because array size is |
|
| //| used to automatically determine FFT length |
|
| //+------------------------------------------------------------------+
|
| void CFastFourierTransform::FFTR1DInv(complex &f[],const int n,double &a[])
|
| {
|
| //--- create a variable
|
| int i=0;
|
| //--- create arrays
|
| double h[];
|
| complex fh[];
|
| //--- check
|
| if(!CAp::Assert(n>0,__FUNCTION__+": incorrect N!"))
|
| return;
|
| //--- check
|
| if(!CAp::Assert(CAp::Len(f)>=(int)MathFloor((double)n/2.0)+1,__FUNCTION__+": Length(F)<Floor(N/2)+1!"))
|
| return;
|
| //--- check
|
| if(!CAp::Assert(CMath::IsFinite(f[0].real),__FUNCTION__+": F contains infinite or NAN values!"))
|
| return;
|
| for(i=1; i<=(int)MathFloor((double)n/2.0)-1; i++)
|
| {
|
| //--- check
|
| if(!CAp::Assert(CMath::IsFinite(f[i].real) && CMath::IsFinite(f[i].imag),__FUNCTION__+": F contains infinite or NAN values!"))
|
| return;
|
| }
|
| //--- check
|
| if(!CAp::Assert(CMath::IsFinite(f[(int)MathFloor((double)n/2.0)].real),__FUNCTION__+": F contains infinite or NAN values!"))
|
| return;
|
| //--- check
|
| if(n%2!=0)
|
| {
|
| //--- check
|
| if(!CAp::Assert(CMath::IsFinite(f[(int)MathFloor((double)n/2.0)].imag),__FUNCTION__+": F contains infinite or NAN values!"))
|
| return;
|
| }
|
| //--- Special case: N=1,FFT is just identity transform.
|
| //--- After this block we assume that N is strictly greater than 1.
|
| if(n==1)
|
| {
|
| //--- allocation
|
| ArrayResize(a,1);
|
| a[0]=f[0].real;
|
| //--- exit the function
|
| return;
|
| }
|
| //--- inverse real FFT is reduced to the inverse real FHT,
|
| //--- which is reduced to the forward real FHT,
|
| //--- which is reduced to the forward real FFT.
|
| //--- Don't worry,it is really compact and efficient reduction :)
|
| ArrayResize(h,n);
|
| ArrayResize(a,n);
|
| //--- change values
|
| h[0]=f[0].real;
|
| for(i=1; i<=(int)MathFloor((double)n/2.0)-1; i++)
|
| {
|
| h[i]=f[i].real-f[i].imag;
|
| h[n-i]=f[i].real+f[i].imag;
|
| }
|
| //--- check
|
| if(n%2==0)
|
| h[(int)MathFloor((double)n/2.0)]=f[(int)MathFloor((double)n/2.0)].real;
|
| else
|
| {
|
| h[(int)MathFloor((double)n/2.0)]=f[(int)MathFloor((double)n/2.0)].real-f[(int)MathFloor((double)n/2.0)].imag;
|
| h[(int)MathFloor((double)n/2.0)+1]=f[(int)MathFloor((double)n/2.0)].real+f[(int)MathFloor((double)n/2.0)].imag;
|
| }
|
| //--- function call
|
| FFTR1D(h,n,fh);
|
| //--- change values
|
| for(i=0; i<=n-1; i++)
|
| a[i]=(fh[i].real-fh[i].imag)/n;
|
| }
|
| //+------------------------------------------------------------------+
|
| //| |
|
| //+------------------------------------------------------------------+
|
| void CFastFourierTransform::FFTR1DInv(CRowComplex &f,const int n,CRowDouble &a)
|
| {
|
| //--- create a variable
|
| int i=0;
|
| //--- create arrays
|
| double h[];
|
| complex fh[];
|
| //--- check
|
| if(!CAp::Assert(n>0,__FUNCTION__+": incorrect N!"))
|
| return;
|
| //--- check
|
| if(!CAp::Assert(CAp::Len(f)>=(int)MathFloor((double)n/2.0)+1,__FUNCTION__+": Length(F)<Floor(N/2)+1!"))
|
| return;
|
| //--- check
|
| if(!CAp::Assert(CMath::IsFinite(f[0].real),__FUNCTION__+": F contains infinite or NAN values!"))
|
| return;
|
| for(i=1; i<=(int)MathFloor((double)n/2.0)-1; i++)
|
| {
|
| //--- check
|
| if(!CAp::Assert(CMath::IsFinite(f[i].real) && CMath::IsFinite(f[i].imag),__FUNCTION__+": F contains infinite or NAN values!"))
|
| return;
|
| }
|
| //--- check
|
| if(!CAp::Assert(CMath::IsFinite(f[(int)MathFloor((double)n/2.0)].real),__FUNCTION__+": F contains infinite or NAN values!"))
|
| return;
|
| //--- check
|
| if(n%2!=0)
|
| {
|
| //--- check
|
| if(!CAp::Assert(CMath::IsFinite(f[(int)MathFloor((double)n/2.0)].imag),__FUNCTION__+": F contains infinite or NAN values!"))
|
| return;
|
| }
|
| //--- Special case: N=1,FFT is just identity transform.
|
| //--- After this block we assume that N is strictly greater than 1.
|
| if(n==1)
|
| {
|
| //--- allocation
|
| a.Resize(1);
|
| a.Set(0,f[0].real);
|
| //--- exit the function
|
| return;
|
| }
|
| //--- inverse real FFT is reduced to the inverse real FHT,
|
| //--- which is reduced to the forward real FHT,
|
| //--- which is reduced to the forward real FFT.
|
| //--- Don't worry,it is really compact and efficient reduction :)
|
| ArrayResize(h,n);
|
| a.Resize(n);
|
| //--- change values
|
| h[0]=f[0].real;
|
| for(i=1; i<=(int)MathFloor((double)n/2.0)-1; i++)
|
| {
|
| h[i]=f[i].real-f[i].imag;
|
| h[n-i]=f[i].real+f[i].imag;
|
| }
|
| //--- check
|
| if(n%2==0)
|
| h[(int)MathFloor((double)n/2.0)]=f[(int)MathFloor((double)n/2.0)].real;
|
| else
|
| {
|
| h[(int)MathFloor((double)n/2.0)]=f[(int)MathFloor((double)n/2.0)].real-f[(int)MathFloor((double)n/2.0)].imag;
|
| h[(int)MathFloor((double)n/2.0)+1]=f[(int)MathFloor((double)n/2.0)].real+f[(int)MathFloor((double)n/2.0)].imag;
|
| }
|
| //--- function call
|
| FFTR1D(h,n,fh);
|
| //--- change values
|
| for(i=0; i<n; i++)
|
| a.Set(i,(fh[i].real-fh[i].imag)/n);
|
| }
|
| //+------------------------------------------------------------------+
|
| //| Internal subroutine. Never call it directly! |
|
| //+------------------------------------------------------------------+
|
| void CFastFourierTransform::FFTR1DInternalEven(double &a[],const int n,
|
| double &buf[],CFtPlan &plan)
|
| {
|
| CRowDouble A=a;
|
| CRowDouble Buf;
|
| FFTR1DInternalEven(A,n,Buf,plan);
|
| A.ToArray(a);
|
| Buf.ToArray(buf);
|
| }
|
| //+------------------------------------------------------------------+
|
| //| |
|
| //+------------------------------------------------------------------+
|
| void CFastFourierTransform::FFTR1DInternalEven(CRowDouble &a,const int n,
|
| CRowDouble &buf,CFtPlan &plan)
|
| {
|
| //--- create variables
|
| double x=0;
|
| double y=0;
|
| int i=0;
|
| int n2=0;
|
| int idx=0;
|
| complex hn=0;
|
| complex hmnc=0;
|
| complex v=0;
|
| int i_=0;
|
| //--- check
|
| if(!CAp::Assert(n>0 && n%2==0,__FUNCTION__+": incorrect N!"))
|
| return;
|
| //--- Special cases:
|
| //--- * N=2
|
| //--- After this block we assume that N is strictly greater than 2
|
| if(n==2)
|
| {
|
| //--- change values
|
| x=a[0]+a[1];
|
| y=a[0]-a[1];
|
| a.Set(0,x);
|
| a.Set(1,y);
|
| //--- exit the function
|
| return;
|
| }
|
| //--- even-size real FFT,use reduction to the complex task
|
| n2=n/2;
|
| //--- copy
|
| buf=a;
|
| buf.Resize(n);
|
| //--- function call
|
| CFtBase::FtApplyPlan(plan,buf,0,1);
|
| //--- change value
|
| a.Set(0,buf[0]+buf[1]);
|
| for(i=1; i<=n2-1; i++)
|
| {
|
| //--- calculation
|
| idx=2*(i%n2);
|
| hn.real=buf[idx+0];
|
| hn.imag=buf[idx+1];
|
| idx=2*(n2-i);
|
| hmnc.real=buf[idx+0];
|
| hmnc.imag=-buf[idx+1];
|
| v.real=-MathSin(-(2*M_PI*i/n));
|
| v.imag=MathCos(-(2*M_PI*i/n));
|
| v=hn+hmnc-v*(hn-hmnc);
|
| a.Set(2*i+0,0.5*v.real);
|
| a.Set(2*i+1,0.5*v.imag);
|
| }
|
| //--- change value
|
| a.Set(1,buf[0]-buf[1]);
|
| }
|
| //+------------------------------------------------------------------+
|
| //| Internal subroutine. Never call it directly! |
|
| //+------------------------------------------------------------------+
|
| void CFastFourierTransform::FFTR1DInvInternalEven(double &a[],const int n,
|
| double &buf[],CFtPlan &plan)
|
| {
|
| //--- create variables
|
| double x=0;
|
| double y=0;
|
| double t=0;
|
| int i=0;
|
| int n2=0;
|
| //--- check
|
| if(!CAp::Assert(n>0 && n%2==0,__FUNCTION__+": incorrect N!"))
|
| return;
|
| //--- Special cases:
|
| //--- * N=2
|
| //--- After this block we assume that N is strictly greater than 2
|
| if(n==2)
|
| {
|
| //--- change values
|
| x=0.5*(a[0]+a[1]);
|
| y=0.5*(a[0]-a[1]);
|
| a[0]=x;
|
| a[1]=y;
|
| //--- exit the function
|
| return;
|
| }
|
| //--- inverse real FFT is reduced to the inverse real FHT,
|
| //--- which is reduced to the forward real FHT,
|
| //--- which is reduced to the forward real FFT.
|
| //--- Don't worry,it is really compact and efficient reduction :)
|
| n2=n/2;
|
| buf[0]=a[0];
|
| //--- calculation
|
| for(i=1; i<=n2-1; i++)
|
| {
|
| x=a[2*i+0];
|
| y=a[2*i+1];
|
| buf[i]=x-y;
|
| buf[n-i]=x+y;
|
| }
|
| buf[n2]=a[1];
|
| //--- function call
|
| FFTR1DInternalEven(buf,n,a,plan);
|
| //--- change values
|
| a[0]=buf[0]/n;
|
| t=1.0/(double)n;
|
| //--- calculation
|
| for(i=1; i<=n2-1; i++)
|
| {
|
| x=buf[2*i+0];
|
| y=buf[2*i+1];
|
| a[i]=t*(x-y);
|
| a[n-i]=t*(x+y);
|
| }
|
| //--- change value
|
| a[n2]=buf[1]/n;
|
| }
|
| //+------------------------------------------------------------------+
|
| //| Internal subroutine. Never call it directly! |
|
| //+------------------------------------------------------------------+
|
| void CFastFourierTransform::FFTR1DInvInternalEven(CRowDouble &a,const int n,
|
| CRowDouble &buf,CFtPlan &plan)
|
| {
|
| //--- create variables
|
| double x=0;
|
| double y=0;
|
| double t=0;
|
| int i=0;
|
| int n2=0;
|
| //--- check
|
| if(!CAp::Assert(n>0 && n%2==0,__FUNCTION__+": incorrect N!"))
|
| return;
|
| //--- Special cases:
|
| //--- * N=2
|
| //--- After this block we assume that N is strictly greater than 2
|
| if(n==2)
|
| {
|
| //--- change values
|
| x=0.5*(a[0]+a[1]);
|
| y=0.5*(a[0]-a[1]);
|
| a.Set(0,x);
|
| a.Set(1,y);
|
| //--- exit the function
|
| return;
|
| }
|
| //--- inverse real FFT is reduced to the inverse real FHT,
|
| //--- which is reduced to the forward real FHT,
|
| //--- which is reduced to the forward real FFT.
|
| //--- Don't worry,it is really compact and efficient reduction :)
|
| n2=n/2;
|
| buf.Set(0,a[0]);
|
| //--- calculation
|
| for(i=1; i<=n2-1; i++)
|
| {
|
| x=a[2*i+0];
|
| y=a[2*i+1];
|
| buf.Set(i,x-y);
|
| buf.Set(n-i,x+y);
|
| }
|
| buf.Set(n2,a[1]);
|
| //--- function call
|
| FFTR1DInternalEven(buf,n,a,plan);
|
| //--- change values
|
| a.Set(0,buf[0]/n);
|
| t=1.0/(double)n;
|
| //--- calculation
|
| for(i=1; i<=n2-1; i++)
|
| {
|
| x=buf[2*i+0];
|
| y=buf[2*i+1];
|
| a.Set(i,t*(x-y));
|
| a.Set(n-i,t*(x+y));
|
| }
|
| //--- change value
|
| a.Set(n2,buf[1]/n);
|
| }
|
| //+------------------------------------------------------------------+
|
| //| Convolution class |
|
| //+------------------------------------------------------------------+
|
| class CConv
|
| {
|
| public:
|
| static void ConvC1D(complex &a[],const int m,complex &b[],const int n,complex &r[]);
|
| static void ConvC1DInv(complex &a[],const int m,complex &b[],const int n,complex &r[]);
|
| static void ConvC1DCircular(complex &s[],const int m,complex &r[],const int n,complex &c[]);
|
| static void ConvC1DCircularInv(complex &a[],const int m,complex &b[],const int n,complex &r[]);
|
| static void ConvR1D(double &a[],const int m,double &b[],const int n,double &r[]);
|
| static void ConvR1DInv(double &a[],const int m,double &b[],const int n,double &r[]);
|
| static void ConvR1DCircular(double &s[],const int m,double &r[],const int n,double &c[]);
|
| static void ConvR1DCircularInv(double &a[],const int m,double &b[],const int n,double &r[]);
|
| static void ConvC1DX(complex &a[],const int m,complex &b[],const int n,const bool circular,int alg,int q,complex &r[]);
|
| static void ConvR1DX(double &a[],const int m,double &b[],const int n,const bool circular,int alg,int q,double &r[]);
|
| };
|
| //+------------------------------------------------------------------+
|
| //| 1-dimensional complex convolution. |
|
| //| For given A/B returns conv(A,B) (non-circular). Subroutine can |
|
| //| automatically choose between three implementations: |
|
| //| straightforward O(M*N) formula for very small N (or M), |
|
| //| significantly larger than min(M,N), but O(M*N) algorithm is too |
|
| //| slow, and general FFT-based formula for cases where two previois |
|
| //| algorithms are too slow. |
|
| //| Algorithm has max(M,N)*log(max(M,N)) complexity for any M/N. |
|
| //| INPUT PARAMETERS |
|
| //| A - array[0..M-1] - complex function to be transformed |
|
| //| M - problem size |
|
| //| B - array[0..N-1] - complex function to be transformed |
|
| //| N - problem size |
|
| //| OUTPUT PARAMETERS |
|
| //| R - convolution: A*B. array[0..N+M-2]. |
|
| //| NOTE: |
|
| //| It is assumed that A is zero at T<0, B is zero too. If one or|
|
| //| both functions have non-zero values at negative T's, you can |
|
| //| still use this subroutine - just shift its result |
|
| //| correspondingly. |
|
| //+------------------------------------------------------------------+
|
| void CConv::ConvC1D(complex &a[],const int m,complex &b[],const int n,complex &r[])
|
| {
|
| //--- check
|
| if(!CAp::Assert(n>0 && m>0,__FUNCTION__+": incorrect N or M!"))
|
| return;
|
| //--- normalize task: make M>=N,
|
| //--- so A will be longer that B.
|
| if(m<n)
|
| {
|
| //--- function call
|
| ConvC1D(b,n,a,m,r);
|
| return;
|
| }
|
| //--- function call
|
| ConvC1DX(a,m,b,n,false,-1,0,r);
|
| }
|
| //+------------------------------------------------------------------+
|
| //| 1-dimensional complex non-circular deconvolution (inverse of |
|
| //| ConvC1D()). |
|
| //| Algorithm has M*log(M)) complexity for any M (composite or prime)|
|
| //| INPUT PARAMETERS |
|
| //| A - array[0..M-1] - convolved signal, A = conv(R, B) |
|
| //| M - convolved signal length |
|
| //| B - array[0..N-1] - response |
|
| //| N - response length, N<=M |
|
| //| OUTPUT PARAMETERS |
|
| //| R - deconvolved signal. array[0..M-N]. |
|
| //| NOTE: |
|
| //| deconvolution is unstable process and may result in division |
|
| //| by zero (if your response function is degenerate, i.e. has zero |
|
| //| Fourier coefficient). |
|
| //| NOTE: |
|
| //| It is assumed that A is zero at T<0, B is zero too. If one |
|
| //| or both functions have non-zero values at negative T's, you can |
|
| //| still use this subroutine - just shift its result correspondingly|
|
| //+------------------------------------------------------------------+
|
| void CConv::ConvC1DInv(complex &a[],const int m,complex &b[],
|
| const int n,complex &r[])
|
| {
|
| //--- create variables
|
| int i=0;
|
| int p=0;
|
| complex c1=0;
|
| complex c2=0;
|
| complex c3=0;
|
| double t=0;
|
| //--- create arrays
|
| double buf[];
|
| double buf2[];
|
| //--- object of class
|
| CFtPlan plan;
|
| //--- check
|
| if(!CAp::Assert((n>0 && m>0) && n<=m,__FUNCTION__+": incorrect N or M!"))
|
| return;
|
| //--- function call
|
| p=CFtBase::FtBaseFindSmooth(m);
|
| CFtBase::FtComplexFFTPlan(p,1,plan);
|
| //--- allocation
|
| ArrayResize(buf,2*p);
|
| //--- copy
|
| for(i=0; i<=m-1; i++)
|
| {
|
| buf[2*i+0]=a[i].real;
|
| buf[2*i+1]=a[i].imag;
|
| }
|
| //--- make zero
|
| for(i=m; i<=p-1; i++)
|
| {
|
| buf[2*i+0]=0;
|
| buf[2*i+1]=0;
|
| }
|
| //--- allocation
|
| ArrayResize(buf2,2*p);
|
| //--- copy
|
| for(i=0; i<=n-1; i++)
|
| {
|
| buf2[2*i+0]=b[i].real;
|
| buf2[2*i+1]=b[i].imag;
|
| }
|
| //--- make zero
|
| for(i=n; i<=p-1; i++)
|
| {
|
| buf2[2*i+0]=0;
|
| buf2[2*i+1]=0;
|
| }
|
| //--- function call
|
| CFtBase::FtApplyPlan(plan,buf,0,1);
|
| CFtBase::FtApplyPlan(plan,buf2,0,1);
|
| //--- calculation
|
| for(i=0; i<=p-1; i++)
|
| {
|
| c1.real=buf[2*i+0];
|
| c1.imag=buf[2*i+1];
|
| c2.real=buf2[2*i+0];
|
| c2.imag=buf2[2*i+1];
|
| c3=c1/c2;
|
| buf[2*i+0]=c3.real;
|
| buf[2*i+1]=-c3.imag;
|
| }
|
| //--- function call
|
| CFtBase::FtApplyPlan(plan,buf,0,1);
|
| t=1.0/(double)p;
|
| //--- allocation
|
| ArrayResize(r,m-n+1);
|
| //--- change values
|
| for(i=0; i<=m-n; i++)
|
| {
|
| r[i].real=t*buf[2*i+0];
|
| r[i].imag=-(t*buf[2*i+1]);
|
| }
|
| }
|
| //+------------------------------------------------------------------+
|
| //| 1-dimensional circular complex convolution. |
|
| //| For given S/R returns conv(S,R) (circular). Algorithm has |
|
| //| linearithmic complexity for any M/N. |
|
| //| IMPORTANT: normal convolution is commutative, i.e. it is |
|
| //| symmetric - conv(A,B)=conv(B,A). Cyclic convolution IS NOT. One |
|
| //| function - S - is a signal, periodic function, and another - R - |
|
| //| is a response, non-periodic function with limited length. |
|
| //| INPUT PARAMETERS |
|
| //| S - array[0..M-1] - complex periodic signal |
|
| //| M - problem size |
|
| //| B - array[0..N-1] - complex non-periodic response |
|
| //| N - problem size |
|
| //| OUTPUT PARAMETERS |
|
| //| R - convolution: A*B. array[0..M-1]. |
|
| //| NOTE: |
|
| //| It is assumed that B is zero at T<0. If it has non-zero |
|
| //| values at negative T's, you can still use this subroutine - just |
|
| //| shift its result correspondingly. |
|
| //+------------------------------------------------------------------+
|
| void CConv::ConvC1DCircular(complex &s[],const int m,complex &r[],
|
| const int n,complex &c[])
|
| {
|
| //--- create variables
|
| int i1=0;
|
| int i2=0;
|
| int j2=0;
|
| int i_=0;
|
| int i1_=0;
|
| //--- create array
|
| complex buf[];
|
| //--- check
|
| if(!CAp::Assert(n>0 && m>0,__FUNCTION__+": incorrect N or M!"))
|
| return;
|
| //--- normalize task: make M>=N,
|
| //--- so A will be longer (at least - not shorter) that B.
|
| if(m<n)
|
| {
|
| //--- allocation
|
| ArrayResize(buf,m);
|
| //--- initialization
|
| for(i1=0; i1<=m-1; i1++)
|
| buf[i1]=0;
|
| i1=0;
|
| //--- calculation
|
| while(i1<n)
|
| {
|
| //--- change values
|
| i2=MathMin(i1+m-1,n-1);
|
| j2=i2-i1;
|
| i1_=i1;
|
| for(i_=0; i_<=j2; i_++)
|
| buf[i_]=buf[i_]+r[i_+i1_];
|
| i1=i1+m;
|
| }
|
| //--- function call
|
| ConvC1DCircular(s,m,buf,m,c);
|
| //--- exit the function
|
| return;
|
| }
|
| //--- function call
|
| ConvC1DX(s,m,r,n,true,-1,0,c);
|
| }
|
| //+------------------------------------------------------------------+
|
| //| 1-dimensional circular complex deconvolution (inverse of |
|
| //| ConvC1DCircular()). |
|
| //| Algorithm has M*log(M)) complexity for any M (composite or prime)|
|
| //| INPUT PARAMETERS |
|
| //| A - array[0..M-1] - convolved periodic signal, |
|
| //| A = conv(R, B) |
|
| //| M - convolved signal length |
|
| //| B - array[0..N-1] - non-periodic response |
|
| //| N - response length |
|
| //| OUTPUT PARAMETERS |
|
| //| R - deconvolved signal. array[0..M-1]. |
|
| //| NOTE: |
|
| //| deconvolution is unstable process and may result in division |
|
| //| by zero (if your response function is degenerate, i.e. has zero |
|
| //| Fourier coefficient). |
|
| //| NOTE: |
|
| //| It is assumed that B is zero at T<0. If it has non-zero |
|
| //| values at negative T's, you can still use this subroutine - just |
|
| //| shift its result correspondingly. |
|
| //+------------------------------------------------------------------+
|
| void CConv::ConvC1DCircularInv(complex &a[],const int m,complex &b[],
|
| const int n,complex &r[])
|
| {
|
| //--- create variables
|
| int i=0;
|
| int i1=0;
|
| int i2=0;
|
| int j2=0;
|
| complex c1=0;
|
| complex c2=0;
|
| complex c3=0;
|
| double t=0;
|
| int i_=0;
|
| int i1_=0;
|
| //--- create arrays
|
| double buf[];
|
| double buf2[];
|
| complex cbuf[];
|
| //--- object of class
|
| CFtPlan plan;
|
| //--- check
|
| if(!CAp::Assert(n>0 && m>0,__FUNCTION__+": incorrect N or M!"))
|
| return;
|
| //--- normalize task: make M>=N,
|
| //--- so A will be longer (at least - not shorter) that B.
|
| if(m<n)
|
| {
|
| //--- allocation
|
| ArrayResize(cbuf,m);
|
| //--- initialization
|
| for(i=0; i<=m-1; i++)
|
| cbuf[i]=0;
|
| i1=0;
|
| //--- calculation
|
| while(i1<n)
|
| {
|
| //--- change values
|
| i2=MathMin(i1+m-1,n-1);
|
| j2=i2-i1;
|
| i1_=i1;
|
| for(i_=0; i_<=j2; i_++)
|
| cbuf[i_]=cbuf[i_]+b[i_+i1_];
|
| i1=i1+m;
|
| }
|
| //--- function call
|
| ConvC1DCircularInv(a,m,cbuf,m,r);
|
| //--- exit the function
|
| return;
|
| }
|
| //--- Task is normalized
|
| CFtBase::FtComplexFFTPlan(m,1,plan);
|
| //--- allocation
|
| ArrayResize(buf,2*m);
|
| //--- copy
|
| for(i=0; i<=m-1; i++)
|
| {
|
| buf[2*i+0]=a[i].real;
|
| buf[2*i+1]=a[i].imag;
|
| }
|
| //--- allocation
|
| ArrayResize(buf2,2*m);
|
| //--- copy
|
| for(i=0; i<=n-1; i++)
|
| {
|
| buf2[2*i+0]=b[i].real;
|
| buf2[2*i+1]=b[i].imag;
|
| }
|
| //--- make zero
|
| for(i=n; i<=m-1; i++)
|
| {
|
| buf2[2*i+0]=0;
|
| buf2[2*i+1]=0;
|
| }
|
| //--- function call
|
| CFtBase::FtApplyPlan(plan,buf,0,1);
|
| CFtBase::FtApplyPlan(plan,buf2,0,1);
|
| //--- calculation
|
| for(i=0; i<=m-1; i++)
|
| {
|
| c1.real=buf[2*i+0];
|
| c1.imag=buf[2*i+1];
|
| c2.real=buf2[2*i+0];
|
| c2.imag=buf2[2*i+1];
|
| c3=c1/c2;
|
| buf[2*i+0]=c3.real;
|
| buf[2*i+1]=-c3.imag;
|
| }
|
| //--- function call
|
| CFtBase::FtApplyPlan(plan,buf,0,1);
|
| t=1.0/(double)m;
|
| //--- allocation
|
| ArrayResize(r,m);
|
| //--- change values
|
| for(i=0; i<=m-1; i++)
|
| {
|
| r[i].real=t*buf[2*i+0];
|
| r[i].imag=-(t*buf[2*i+1]);
|
| }
|
| }
|
| //+------------------------------------------------------------------+
|
| //| 1-dimensional real convolution. |
|
| //| Analogous to ConvC1D(), see ConvC1D() comments for more details. |
|
| //| INPUT PARAMETERS |
|
| //| A - array[0..M-1] - real function to be transformed |
|
| //| M - problem size |
|
| //| B - array[0..N-1] - real function to be transformed |
|
| //| N - problem size |
|
| //| OUTPUT PARAMETERS |
|
| //| R - convolution: A*B. array[0..N+M-2]. |
|
| //| NOTE: |
|
| //| It is assumed that A is zero at T<0, B is zero too. If one |
|
| //| or both functions have non-zero values at negative T's, you can |
|
| //| still use this subroutine - just shift its result correspondingly|
|
| //+------------------------------------------------------------------+
|
| void CConv::ConvR1D(double &a[],const int m,double &b[],
|
| const int n,double &r[])
|
| {
|
| //--- check
|
| if(!CAp::Assert(n>0 && m>0,__FUNCTION__+": incorrect N or M!"))
|
| return;
|
| //--- normalize task: make M>=N,
|
| //--- so A will be longer that B.
|
| if(m<n)
|
| {
|
| //--- function call
|
| ConvR1D(b,n,a,m,r);
|
| return;
|
| }
|
| //--- function call
|
| ConvR1DX(a,m,b,n,false,-1,0,r);
|
| }
|
| //+------------------------------------------------------------------+
|
| //| 1-dimensional real deconvolution (inverse of ConvC1D()). |
|
| //| Algorithm has M*log(M)) complexity for any M (composite or prime)|
|
| //| INPUT PARAMETERS |
|
| //| A - array[0..M-1] - convolved signal, A = conv(R, B) |
|
| //| M - convolved signal length |
|
| //| B - array[0..N-1] - response |
|
| //| N - response length, N<=M |
|
| //| OUTPUT PARAMETERS |
|
| //| R - deconvolved signal. array[0..M-N]. |
|
| //| NOTE: |
|
| //| deconvolution is unstable process and may result in division |
|
| //| by zero (if your response function is degenerate, i.e. has zero |
|
| //| Fourier coefficient). |
|
| //| NOTE: |
|
| //| It is assumed that A is zero at T<0, B is zero too. If one or|
|
| //| both functions have non-zero values at negative T's, you can |
|
| //| still use this subroutine - just shift its result correspondingly|
|
| //+------------------------------------------------------------------+
|
| void CConv::ConvR1DInv(double &a[],const int m,double &b[],
|
| const int n,double &r[])
|
| {
|
| //--- create variables
|
| int i=0;
|
| int p=0;
|
| complex c1=0;
|
| complex c2=0;
|
| complex c3=0;
|
| int i_=0;
|
| //--- create arrays
|
| double buf[];
|
| double buf2[];
|
| double buf3[];
|
| //--- object of class
|
| CFtPlan plan;
|
| //--- check
|
| if(!CAp::Assert((n>0 && m>0) && n<=m,__FUNCTION__+": incorrect N or M!"))
|
| return;
|
| //--- function call
|
| p=CFtBase::FtBaseFindSmoothEven(m);
|
| //--- allocation
|
| ArrayResize(buf,p);
|
| //--- copy
|
| for(i_=0; i_<=m-1; i_++)
|
| buf[i_]=a[i_];
|
| //--- make zero
|
| for(i=m; i<=p-1; i++)
|
| buf[i]=0;
|
| //--- allocation
|
| ArrayResize(buf2,p);
|
| for(i_=0; i_<=n-1; i_++)
|
| buf2[i_]=b[i_];
|
| //--- make zero
|
| for(i=n; i<=p-1; i++)
|
| buf2[i]=0;
|
| //--- allocation
|
| ArrayResize(buf3,p);
|
| //--- function call
|
| CFtBase::FtComplexFFTPlan(p/2,1,plan);
|
| //--- function call
|
| CFastFourierTransform::FFTR1DInternalEven(buf,p,buf3,plan);
|
| //--- function call
|
| CFastFourierTransform::FFTR1DInternalEven(buf2,p,buf3,plan);
|
| //--- change values
|
| buf[0]=buf[0]/buf2[0];
|
| buf[1]=buf[1]/buf2[1];
|
| //--- calculation
|
| for(i=1; i<=p/2-1; i++)
|
| {
|
| c1.real=buf[2*i+0];
|
| c1.imag=buf[2*i+1];
|
| c2.real=buf2[2*i+0];
|
| c2.imag=buf2[2*i+1];
|
| c3=c1/c2;
|
| buf[2*i+0]=c3.real;
|
| buf[2*i+1]=c3.imag;
|
| }
|
| //--- function call
|
| CFastFourierTransform::FFTR1DInvInternalEven(buf,p,buf3,plan);
|
| //--- allocation
|
| ArrayResize(r,m-n+1);
|
| //--- copy
|
| for(i_=0; i_<=m-n; i_++)
|
| r[i_]=buf[i_];
|
| }
|
| //+------------------------------------------------------------------+
|
| //| 1-dimensional circular real convolution. |
|
| //| Analogous to ConvC1DCircular(), see ConvC1DCircular() comments |
|
| //| for more details. |
|
| //| INPUT PARAMETERS |
|
| //| S - array[0..M-1] - real signal |
|
| //| M - problem size |
|
| //| B - array[0..N-1] - real response |
|
| //| N - problem size |
|
| //| OUTPUT PARAMETERS |
|
| //| R - convolution: A*B. array[0..M-1]. |
|
| //| NOTE: |
|
| //| It is assumed that B is zero at T<0. If it has non-zero |
|
| //| values at negative T's, you can still use this subroutine - just |
|
| //| shift its result correspondingly. |
|
| //+------------------------------------------------------------------+
|
| void CConv::ConvR1DCircular(double &s[],const int m,double &r[],
|
| const int n,double &c[])
|
| {
|
| //--- create variables
|
| int i1=0;
|
| int i2=0;
|
| int j2=0;
|
| int i_=0;
|
| int i1_=0;
|
| //--- create array
|
| double buf[];
|
| //--- check
|
| if(!CAp::Assert(n>0 && m>0,__FUNCTION__+": incorrect N or M!"))
|
| return;
|
| //--- normalize task: make M>=N,
|
| //--- so A will be longer (at least - not shorter) that B.
|
| if(m<n)
|
| {
|
| //--- allocation
|
| ArrayResize(buf,m);
|
| //--- initialization
|
| for(i1=0; i1<=m-1; i1++)
|
| buf[i1]=0;
|
| i1=0;
|
| //--- calculation
|
| while(i1<n)
|
| {
|
| i2=MathMin(i1+m-1,n-1);
|
| j2=i2-i1;
|
| i1_=i1;
|
| for(i_=0; i_<=j2; i_++)
|
| buf[i_]=buf[i_]+r[i_+i1_];
|
| i1=i1+m;
|
| }
|
| ConvR1DCircular(s,m,buf,m,c);
|
| //--- exit the function
|
| return;
|
| }
|
| //--- reduce to usual convolution
|
| ConvR1DX(s,m,r,n,true,-1,0,c);
|
| }
|
| //+------------------------------------------------------------------+
|
| //| 1-dimensional complex deconvolution (inverse of ConvC1D()). |
|
| //| Algorithm has M*log(M)) complexity for any M (composite or prime)|
|
| //| INPUT PARAMETERS |
|
| //| A - array[0..M-1] - convolved signal, A = conv(R, B) |
|
| //| M - convolved signal length |
|
| //| B - array[0..N-1] - response |
|
| //| N - response length |
|
| //| OUTPUT PARAMETERS |
|
| //| R - deconvolved signal. array[0..M-N]. |
|
| //| NOTE: |
|
| //| deconvolution is unstable process and may result in division |
|
| //| by zero (if your response function is degenerate, i.e. has zero |
|
| //| Fourier coefficient). |
|
| //| NOTE: |
|
| //| It is assumed that B is zero at T<0. If it has non-zero |
|
| //| values at negative T's, you can still use this subroutine - just |
|
| //| shift its result correspondingly. |
|
| //+------------------------------------------------------------------+
|
| void CConv::ConvR1DCircularInv(double &a[],const int m,double &b[],
|
| const int n,double &r[])
|
| {
|
| //--- create variables
|
| int i=0;
|
| int i1=0;
|
| int i2=0;
|
| int j2=0;
|
| complex c1=0;
|
| complex c2=0;
|
| complex c3=0;
|
| int i_=0;
|
| int i1_=0;
|
| //--- create arrays
|
| double buf[];
|
| double buf2[];
|
| double buf3[];
|
| complex cbuf[];
|
| complex cbuf2[];
|
| //--- object of class
|
| CFtPlan plan;
|
| //--- check
|
| if(!CAp::Assert(n>0 && m>0,__FUNCTION__+": incorrect N or M!"))
|
| return;
|
| //--- normalize task: make M>=N,
|
| //--- so A will be longer (at least - not shorter) that B.
|
| if(m<n)
|
| {
|
| //--- allocation
|
| ArrayResize(buf,m);
|
| //--- initialization
|
| for(i=0; i<=m-1; i++)
|
| buf[i]=0;
|
| i1=0;
|
| //--- calculation
|
| while(i1<n)
|
| {
|
| i2=MathMin(i1+m-1,n-1);
|
| j2=i2-i1;
|
| i1_=i1;
|
| for(i_=0; i_<=j2; i_++)
|
| buf[i_]=buf[i_]+b[i_+i1_];
|
| i1=i1+m;
|
| }
|
| //--- function call
|
| ConvR1DCircularInv(a,m,buf,m,r);
|
| //--- exit the function
|
| return;
|
| }
|
| //--- Task is normalized
|
| if(m%2==0)
|
| {
|
| //--- size is even,use fast even-size FFT
|
| ArrayResize(buf,m);
|
| for(i_=0; i_<=m-1; i_++)
|
| buf[i_]=a[i_];
|
| //--- allocation
|
| ArrayResize(buf2,m);
|
| //--- copy
|
| for(i_=0; i_<=n-1; i_++)
|
| buf2[i_]=b[i_];
|
| //--- make zero
|
| for(i=n; i<=m-1; i++)
|
| buf2[i]=0;
|
| //--- allocation
|
| ArrayResize(buf3,m);
|
| //--- function call
|
| CFtBase::FtComplexFFTPlan(m/2,1,plan);
|
| //--- function call
|
| CFastFourierTransform::FFTR1DInternalEven(buf,m,buf3,plan);
|
| //--- function call
|
| CFastFourierTransform::FFTR1DInternalEven(buf2,m,buf3,plan);
|
| //--- change values
|
| buf[0]=buf[0]/buf2[0];
|
| buf[1]=buf[1]/buf2[1];
|
| //--- calculation
|
| for(i=1; i<=m/2-1; i++)
|
| {
|
| c1.real=buf[2*i+0];
|
| c1.imag=buf[2*i+1];
|
| c2.real=buf2[2*i+0];
|
| c2.imag=buf2[2*i+1];
|
| c3=c1/c2;
|
| buf[2*i+0]=c3.real;
|
| buf[2*i+1]=c3.imag;
|
| }
|
| //--- function call
|
| CFastFourierTransform::FFTR1DInvInternalEven(buf,m,buf3,plan);
|
| //--- allocation
|
| ArrayResize(r,m);
|
| //--- copy
|
| for(i_=0; i_<=m-1; i_++)
|
| r[i_]=buf[i_];
|
| }
|
| else
|
| {
|
| //--- odd-size,use general real FFT
|
| CFastFourierTransform::FFTR1D(a,m,cbuf);
|
| //--- allocation
|
| ArrayResize(buf2,m);
|
| //--- copy
|
| for(i_=0; i_<=n-1; i_++)
|
| buf2[i_]=b[i_];
|
| //--- initialization
|
| for(i=n; i<=m-1; i++)
|
| buf2[i]=0;
|
| //--- function call
|
| CFastFourierTransform::FFTR1D(buf2,m,cbuf2);
|
| //--- calculation
|
| for(i=0; i<=(int)MathFloor((double)m/2.0); i++)
|
| cbuf[i]=cbuf[i]/cbuf2[i];
|
| //--- function call
|
| CFastFourierTransform::FFTR1DInv(cbuf,m,r);
|
| }
|
| }
|
| //+------------------------------------------------------------------+
|
| //| 1-dimensional complex convolution. |
|
| //| Extended subroutine which allows to choose convolution algorithm.|
|
| //| Intended for internal use, ALGLIB users should call |
|
| //| ConvC1D()/ConvC1DCircular(). |
|
| //| INPUT PARAMETERS |
|
| //| A - array[0..M-1] - complex function to be transformed |
|
| //| M - problem size |
|
| //| B - array[0..N-1] - complex function to be transformed |
|
| //| N - problem size, N<=M |
|
| //| Alg - algorithm type: |
|
| //| *-2 auto-select Q for overlap-add |
|
| //| *-1 auto-select algorithm and parameters |
|
| //| * 0 straightforward formula for small N's |
|
| //| * 1 general FFT-based code |
|
| //| * 2 overlap-add with length Q |
|
| //| Q - length for overlap-add |
|
| //| OUTPUT PARAMETERS |
|
| //| R - convolution: A*B. array[0..N+M-1]. |
|
| //+------------------------------------------------------------------+
|
| void CConv::ConvC1DX(complex &a[],const int m,complex &b[],const int n,
|
| const bool circular,int alg,int q,complex &r[])
|
| {
|
| //--- create variables
|
| int i=0;
|
| int j=0;
|
| int p=0;
|
| int ptotal=0;
|
| int i1=0;
|
| int i2=0;
|
| int j1=0;
|
| int j2=0;
|
| complex v=0;
|
| double ax=0;
|
| double ay=0;
|
| double bx=0;
|
| double by=0;
|
| double t=0;
|
| double tx=0;
|
| double ty=0;
|
| double flopcand=0;
|
| double flopbest=0;
|
| int algbest=0;
|
| int i_=0;
|
| int i1_=0;
|
| //--- create arrays
|
| complex bbuf[];
|
| double buf[];
|
| double buf2[];
|
| //--- object of class
|
| CFtPlan plan;
|
| //--- check
|
| if(!CAp::Assert(n>0 && m>0,__FUNCTION__+": incorrect N or M!"))
|
| return;
|
| //--- check
|
| if(!CAp::Assert(n<=m,__FUNCTION__+": N<M assumption is false!"))
|
| return;
|
| //--- Auto-select
|
| if(alg==-1 || alg==-2)
|
| {
|
| //--- Initial candidate: straightforward implementation.
|
| //--- If we want to use auto-fitted overlap-add,
|
| //--- flop count is initialized by large real number - to force
|
| //--- another algorithm selection
|
| algbest=0;
|
| //--- check
|
| if(alg==-1)
|
| flopbest=2*m*n;
|
| else
|
| flopbest=CMath::m_maxrealnumber;
|
| //--- Another candidate - generic FFT code
|
| if(alg==-1)
|
| {
|
| //--- check
|
| if(circular && CFtBase::FtBaseIsSmooth(m))
|
| {
|
| //--- special code for circular convolution of a sequence with a smooth length
|
| flopcand=3*CFtBase::FtBaseGetFlopEstimate(m)+6*m;
|
| //--- check
|
| if(flopcand<flopbest)
|
| {
|
| algbest=1;
|
| flopbest=flopcand;
|
| }
|
| }
|
| else
|
| {
|
| //--- general cyclic/non-cyclic convolution
|
| p=CFtBase::FtBaseFindSmooth(m+n-1);
|
| flopcand=3*CFtBase::FtBaseGetFlopEstimate(p)+6*p;
|
| //--- check
|
| if(flopcand<flopbest)
|
| {
|
| algbest=1;
|
| flopbest=flopcand;
|
| }
|
| }
|
| }
|
| //--- Another candidate - overlap-add
|
| q=1;
|
| ptotal=1;
|
| while(ptotal<n)
|
| ptotal=ptotal*2;
|
| //--- calculation
|
| while(ptotal<=m+n-1)
|
| {
|
| //--- change values
|
| p=ptotal-n+1;
|
| flopcand=(int)MathCeil((double)m/(double)p)*(2*CFtBase::FtBaseGetFlopEstimate(ptotal)+8*ptotal);
|
| //--- check
|
| if(flopcand<flopbest)
|
| {
|
| flopbest=flopcand;
|
| algbest=2;
|
| q=p;
|
| }
|
| ptotal=ptotal*2;
|
| }
|
| //--- change value
|
| alg=algbest;
|
| //--- function call
|
| ConvC1DX(a,m,b,n,circular,alg,q,r);
|
| //--- exit the function
|
| return;
|
| }
|
| //--- straightforward formula for
|
| //--- circular and non-circular convolutions.
|
| //--- Very simple code,no further comments needed.
|
| if(alg==0)
|
| {
|
| //--- Special case: N=1
|
| if(n==1)
|
| {
|
| //--- allocation
|
| ArrayResize(r,m);
|
| v=b[0];
|
| for(i_=0; i_<=m-1; i_++)
|
| r[i_]=v*a[i_];
|
| //--- exit the function
|
| return;
|
| }
|
| //--- use straightforward formula
|
| if(circular)
|
| {
|
| //--- circular convolution
|
| ArrayResize(r,m);
|
| //--- change values
|
| v=b[0];
|
| for(i_=0; i_<=m-1; i_++)
|
| r[i_]=v*a[i_];
|
| //--- calculation
|
| for(i=1; i<=n-1; i++)
|
| {
|
| //--- change values
|
| v=b[i];
|
| i1=0;
|
| i2=i-1;
|
| j1=m-i;
|
| j2=m-1;
|
| i1_=j1-i1;
|
| //--- calculation
|
| for(i_=i1; i_<=i2; i_++)
|
| r[i_]=r[i_]+v*a[i_+i1_];
|
| //--- change values
|
| i1=i;
|
| i2=m-1;
|
| j1=0;
|
| j2=m-i-1;
|
| i1_=j1-i1;
|
| //--- calculation
|
| for(i_=i1; i_<=i2; i_++)
|
| r[i_]=r[i_]+v*a[i_+i1_];
|
| }
|
| }
|
| else
|
| {
|
| //--- non-circular convolution
|
| ArrayResize(r,m+n-1);
|
| //--- initialization
|
| for(i=0; i<=m+n-2; i++)
|
| r[i]=0;
|
| //--- calculation
|
| for(i=0; i<=n-1; i++)
|
| {
|
| v=b[i];
|
| i1_=-i;
|
| for(i_=i; i_<=i+m-1; i_++)
|
| r[i_]=r[i_]+v*a[i_+i1_];
|
| }
|
| }
|
| //--- exit the function
|
| return;
|
| }
|
| //--- general FFT-based code for
|
| //--- circular and non-circular convolutions.
|
| //--- First,if convolution is circular,we test whether M is smooth or not.
|
| //--- If it is smooth,we just use M-length FFT to calculate convolution.
|
| //--- If it is not,we calculate non-circular convolution and wrap it arount.
|
| //--- IF convolution is non-circular,we use zero-padding + FFT.
|
| if(alg==1)
|
| {
|
| //--- check
|
| if(circular && CFtBase::FtBaseIsSmooth(m))
|
| {
|
| //--- special code for circular convolution with smooth M
|
| CFtBase::FtComplexFFTPlan(m,1,plan);
|
| //--- allocation
|
| ArrayResize(buf,2*m);
|
| //--- copy
|
| for(i=0; i<=m-1; i++)
|
| {
|
| buf[2*i+0]=a[i].real;
|
| buf[2*i+1]=a[i].imag;
|
| }
|
| //--- allocation
|
| ArrayResize(buf2,2*m);
|
| //--- copy
|
| for(i=0; i<=n-1; i++)
|
| {
|
| buf2[2*i+0]=b[i].real;
|
| buf2[2*i+1]=b[i].imag;
|
| }
|
| //--- make zero
|
| for(i=n; i<=m-1; i++)
|
| {
|
| buf2[2*i+0]=0;
|
| buf2[2*i+1]=0;
|
| }
|
| //--- function call
|
| CFtBase::FtApplyPlan(plan,buf,0,1);
|
| CFtBase::FtApplyPlan(plan,buf2,0,1);
|
| //--- calculation
|
| for(i=0; i<=m-1; i++)
|
| {
|
| ax=buf[2*i+0];
|
| ay=buf[2*i+1];
|
| bx=buf2[2*i+0];
|
| by=buf2[2*i+1];
|
| tx=ax*bx-ay*by;
|
| ty=ax*by+ay*bx;
|
| buf[2*i+0]=tx;
|
| buf[2*i+1]=-ty;
|
| }
|
| //--- function call
|
| CFtBase::FtApplyPlan(plan,buf,0,1);
|
| t=1.0/(double)m;
|
| //--- allocation
|
| ArrayResize(r,m);
|
| //--- change values
|
| for(i=0; i<=m-1; i++)
|
| {
|
| r[i].real=t*buf[2*i+0];
|
| r[i].imag=-(t*buf[2*i+1]);
|
| }
|
| }
|
| else
|
| {
|
| //--- M is non-smooth,general code (circular/non-circular):
|
| //--- * first part is the same for circular and non-circular
|
| //--- convolutions. zero padding,FFTs,inverse FFTs
|
| //--- * second part differs:
|
| //--- * for non-circular convolution we just copy array
|
| //--- * for circular convolution we add array tail to its head
|
| p=CFtBase::FtBaseFindSmooth(m+n-1);
|
| CFtBase::FtComplexFFTPlan(p,1,plan);
|
| //--- allocation
|
| ArrayResize(buf,2*p);
|
| //--- copy
|
| for(i=0; i<=m-1; i++)
|
| {
|
| buf[2*i+0]=a[i].real;
|
| buf[2*i+1]=a[i].imag;
|
| }
|
| //--- make zero
|
| for(i=m; i<=p-1; i++)
|
| {
|
| buf[2*i+0]=0;
|
| buf[2*i+1]=0;
|
| }
|
| //--- allocation
|
| ArrayResize(buf2,2*p);
|
| //--- copy
|
| for(i=0; i<=n-1; i++)
|
| {
|
| buf2[2*i+0]=b[i].real;
|
| buf2[2*i+1]=b[i].imag;
|
| }
|
| //--- make zero
|
| for(i=n; i<=p-1; i++)
|
| {
|
| buf2[2*i+0]=0;
|
| buf2[2*i+1]=0;
|
| }
|
| //--- function call
|
| CFtBase::FtApplyPlan(plan,buf,0,1);
|
| CFtBase::FtApplyPlan(plan,buf2,0,1);
|
| //--- calculation
|
| for(i=0; i<=p-1; i++)
|
| {
|
| ax=buf[2*i+0];
|
| ay=buf[2*i+1];
|
| bx=buf2[2*i+0];
|
| by=buf2[2*i+1];
|
| tx=ax*bx-ay*by;
|
| ty=ax*by+ay*bx;
|
| buf[2*i+0]=tx;
|
| buf[2*i+1]=-ty;
|
| }
|
| //--- function call
|
| CFtBase::FtApplyPlan(plan,buf,0,1);
|
| t=1.0/(double)p;
|
| //--- check
|
| if(circular)
|
| {
|
| //--- circular,add tail to head
|
| ArrayResize(r,m);
|
| //--- copy
|
| for(i=0; i<=m-1; i++)
|
| {
|
| r[i].real=t*buf[2*i+0];
|
| r[i].imag=-(t*buf[2*i+1]);
|
| }
|
| //--- change values
|
| for(i=m; i<=m+n-2; i++)
|
| {
|
| r[i-m].real=r[i-m].real+t*buf[2*i+0];
|
| r[i-m].imag=r[i-m].imag-t*buf[2*i+1];
|
| }
|
| }
|
| else
|
| {
|
| //--- non-circular,just copy
|
| ArrayResize(r,m+n-1);
|
| for(i=0; i<=m+n-2; i++)
|
| {
|
| r[i].real=t*buf[2*i+0];
|
| r[i].imag=-(t*buf[2*i+1]);
|
| }
|
| }
|
| }
|
| //--- exit the function
|
| return;
|
| }
|
| //--- overlap-add method for
|
| //--- circular and non-circular convolutions.
|
| //--- First part of code (separate FFTs of input blocks) is the same
|
| //--- for all types of convolution. Second part (overlapping outputs)
|
| //--- differs for different types of convolution. We just copy output
|
| //--- when convolution is non-circular. We wrap it around,if it is
|
| //--- circular.
|
| if(alg==2)
|
| {
|
| //--- allocation
|
| ArrayResize(buf,2*(q+n-1));
|
| //--- prepare R
|
| if(circular)
|
| {
|
| //--- allocation
|
| ArrayResize(r,m);
|
| for(i=0; i<=m-1; i++)
|
| r[i]=0;
|
| }
|
| else
|
| {
|
| //--- allocation
|
| ArrayResize(r,m+n-1);
|
| for(i=0; i<=m+n-2; i++)
|
| r[i]=0;
|
| }
|
| //--- pre-calculated FFT(B)
|
| ArrayResize(bbuf,q+n-1);
|
| for(i_=0; i_<=n-1; i_++)
|
| bbuf[i_]=b[i_];
|
| for(j=n; j<=q+n-2; j++)
|
| bbuf[j]=0;
|
| //--- function call
|
| CFastFourierTransform::FFTC1D(bbuf,q+n-1);
|
| //--- prepare FFT plan for chunks of A
|
| CFtBase::FtComplexFFTPlan(q+n-1,1,plan);
|
| //--- main overlap-add cycle
|
| i=0;
|
| //--- calculation
|
| while(i<=m-1)
|
| {
|
| p=MathMin(q,m-i);
|
| //--- copy
|
| for(j=0; j<=p-1; j++)
|
| {
|
| buf[2*j+0]=a[i+j].real;
|
| buf[2*j+1]=a[i+j].imag;
|
| }
|
| //--- make zero
|
| for(j=p; j<=q+n-2; j++)
|
| {
|
| buf[2*j+0]=0;
|
| buf[2*j+1]=0;
|
| }
|
| //--- function call
|
| CFtBase::FtApplyPlan(plan,buf,0,1);
|
| //--- calculation
|
| for(j=0; j<=q+n-2; j++)
|
| {
|
| ax=buf[2*j+0];
|
| ay=buf[2*j+1];
|
| bx=bbuf[j].real;
|
| by=bbuf[j].imag;
|
| tx=ax*bx-ay*by;
|
| ty=ax*by+ay*bx;
|
| buf[2*j+0]=tx;
|
| buf[2*j+1]=-ty;
|
| }
|
| //--- function call
|
| CFtBase::FtApplyPlan(plan,buf,0,1);
|
| t=1.0/(double)(q+n-1);
|
| //--- check
|
| if(circular)
|
| {
|
| j1=MathMin(i+p+n-2,m-1)-i;
|
| j2=j1+1;
|
| }
|
| else
|
| {
|
| j1=p+n-2;
|
| j2=j1+1;
|
| }
|
| //--- change values
|
| for(j=0; j<=j1; j++)
|
| {
|
| r[i+j].real=r[i+j].real+buf[2*j+0]*t;
|
| r[i+j].imag=r[i+j].imag-buf[2*j+1]*t;
|
| }
|
| //--- change values
|
| for(j=j2; j<=p+n-2; j++)
|
| {
|
| r[j-j2].real=r[j-j2].real+buf[2*j+0]*t;
|
| r[j-j2].imag=r[j-j2].imag-buf[2*j+1]*t;
|
| }
|
| i=i+p;
|
| }
|
| //--- exit the function
|
| return;
|
| }
|
| }
|
| //+------------------------------------------------------------------+
|
| //| 1-dimensional real convolution. |
|
| //| Extended subroutine which allows to choose convolution algorithm.|
|
| //| Intended for internal use, ALGLIB users should call ConvR1D(). |
|
| //| INPUT PARAMETERS |
|
| //| A - array[0..M-1] - complex function to be transformed |
|
| //| M - problem size |
|
| //| B - array[0..N-1] - complex function to be transformed |
|
| //| N - problem size, N<=M |
|
| //| Alg - algorithm type: |
|
| //| *-2 auto-select Q for overlap-add |
|
| //| *-1 auto-select algorithm and parameters |
|
| //| * 0 straightforward formula for small N's |
|
| //| * 1 general FFT-based code |
|
| //| * 2 overlap-add with length Q |
|
| //| Q - length for overlap-add |
|
| //| OUTPUT PARAMETERS |
|
| //| R - convolution: A*B. array[0..N+M-1]. |
|
| //+------------------------------------------------------------------+
|
| void CConv::ConvR1DX(double &a[],const int m,double &b[],const int n,
|
| const bool circular,int alg,int q,double &r[])
|
| {
|
| //--- create variables
|
| double v=0;
|
| int i=0;
|
| int j=0;
|
| int p=0;
|
| int ptotal=0;
|
| int i1=0;
|
| int i2=0;
|
| int j1=0;
|
| int j2=0;
|
| double ax=0;
|
| double ay=0;
|
| double bx=0;
|
| double by=0;
|
| double tx=0;
|
| double ty=0;
|
| double flopcand=0;
|
| double flopbest=0;
|
| int algbest=0;
|
| int i_=0;
|
| int i1_=0;
|
| //--- create arrays
|
| double buf[];
|
| double buf2[];
|
| double buf3[];
|
| //--- object of class
|
| CFtPlan plan;
|
| //--- check
|
| if(!CAp::Assert(n>0 && m>0,__FUNCTION__+": incorrect N or M!"))
|
| return;
|
| //--- check
|
| if(!CAp::Assert(n<=m,__FUNCTION__+": N<M assumption is false!"))
|
| return;
|
| //--- handle special cases
|
| if(MathMin(m,n)<=2)
|
| alg=0;
|
| //--- Auto-select
|
| if(alg<0)
|
| {
|
| //--- Initial candidate: straightforward implementation.
|
| //--- If we want to use auto-fitted overlap-add,
|
| //--- flop count is initialized by large real number - to force
|
| //--- another algorithm selection
|
| algbest=0;
|
| //--- check
|
| if(alg==-1)
|
| flopbest=0.15*m*n;
|
| else
|
| flopbest=CMath::m_maxrealnumber;
|
| //--- Another candidate - generic FFT code
|
| if(alg==-1)
|
| {
|
| //--- check
|
| if((circular && CFtBase::FtBaseIsSmooth(m)) && m%2==0)
|
| {
|
| //--- special code for circular convolution of a sequence with a smooth length
|
| flopcand=3*CFtBase::FtBaseGetFlopEstimate(m/2)+(double)(6*m)/2.0;
|
| //--- check
|
| if(flopcand<flopbest)
|
| {
|
| algbest=1;
|
| flopbest=flopcand;
|
| }
|
| }
|
| else
|
| {
|
| //--- general cyclic/non-cyclic convolution
|
| p=CFtBase::FtBaseFindSmoothEven(m+n-1);
|
| flopcand=3*CFtBase::FtBaseGetFlopEstimate(p/2)+(double)(6*p)/2.0;
|
| //--- check
|
| if(flopcand<flopbest)
|
| {
|
| algbest=1;
|
| flopbest=flopcand;
|
| }
|
| }
|
| }
|
| //--- Another candidate - overlap-add
|
| q=1;
|
| ptotal=1;
|
| while(ptotal<n)
|
| ptotal=ptotal*2;
|
| //--- calculation
|
| while(ptotal<=m+n-1)
|
| {
|
| //--- change values
|
| p=ptotal-n+1;
|
| flopcand=(int)MathCeil((double)m/(double)p)*(2*CFtBase::FtBaseGetFlopEstimate(ptotal/2)+1*(ptotal/2));
|
| //--- check
|
| if(flopcand<flopbest)
|
| {
|
| flopbest=flopcand;
|
| algbest=2;
|
| q=p;
|
| }
|
| ptotal=ptotal*2;
|
| }
|
| alg=algbest;
|
| //--- function call
|
| ConvR1DX(a,m,b,n,circular,alg,q,r);
|
| //--- exit the function
|
| return;
|
| }
|
| //--- straightforward formula for
|
| //--- circular and non-circular convolutions.
|
| //--- Very simple code,no further comments needed.
|
| if(alg==0)
|
| {
|
| //--- Special case: N=1
|
| if(n==1)
|
| {
|
| //--- allocation
|
| ArrayResize(r,m);
|
| v=b[0];
|
| for(i_=0; i_<=m-1; i_++)
|
| r[i_]=v*a[i_];
|
| //--- exit the function
|
| return;
|
| }
|
| //--- use straightforward formula
|
| if(circular)
|
| {
|
| //--- circular convolution
|
| ArrayResize(r,m);
|
| v=b[0];
|
| for(i_=0; i_<=m-1; i_++)
|
| r[i_]=v*a[i_];
|
| //--- calculation
|
| for(i=1; i<=n-1; i++)
|
| {
|
| //--- change values
|
| v=b[i];
|
| i1=0;
|
| i2=i-1;
|
| j1=m-i;
|
| j2=m-1;
|
| i1_=j1-i1;
|
| //--- calculation
|
| for(i_=i1; i_<=i2; i_++)
|
| r[i_]=r[i_]+v*a[i_+i1_];
|
| //--- change values
|
| i1=i;
|
| i2=m-1;
|
| j1=0;
|
| j2=m-i-1;
|
| i1_=j1-i1;
|
| //--- calculation
|
| for(i_=i1; i_<=i2; i_++)
|
| r[i_]=r[i_]+v*a[i_+i1_];
|
| }
|
| }
|
| else
|
| {
|
| //--- non-circular convolution
|
| ArrayResize(r,m+n-1);
|
| for(i=0; i<=m+n-2; i++)
|
| r[i]=0;
|
| //--- calculation
|
| for(i=0; i<=n-1; i++)
|
| {
|
| v=b[i];
|
| i1_=-i;
|
| for(i_=i; i_<=i+m-1; i_++)
|
| r[i_]=r[i_]+v*a[i_+i1_];
|
| }
|
| }
|
| //--- exit the function
|
| return;
|
| }
|
| //--- general FFT-based code for
|
| //--- circular and non-circular convolutions.
|
| //--- First,if convolution is circular,we test whether M is smooth or not.
|
| //--- If it is smooth,we just use M-length FFT to calculate convolution.
|
| //--- If it is not,we calculate non-circular convolution and wrap it arount.
|
| //--- If convolution is non-circular,we use zero-padding + FFT.
|
| //--- We assume that M+N-1>2 - we should call small case code otherwise
|
| if(alg==1)
|
| {
|
| //--- check
|
| if(!CAp::Assert(m+n-1>2,__FUNCTION__+": internal error!"))
|
| return;
|
| //--- check
|
| if((circular && CFtBase::FtBaseIsSmooth(m)) && m%2==0)
|
| {
|
| //--- special code for circular convolution with smooth even M
|
| ArrayResize(buf,m);
|
| for(i_=0; i_<=m-1; i_++)
|
| buf[i_]=a[i_];
|
| //--- allocation
|
| ArrayResize(buf2,m);
|
| for(i_=0; i_<=n-1; i_++)
|
| buf2[i_]=b[i_];
|
| for(i=n; i<=m-1; i++)
|
| buf2[i]=0;
|
| //--- allocation
|
| ArrayResize(buf3,m);
|
| //--- function call
|
| CFtBase::FtComplexFFTPlan(m/2,1,plan);
|
| //--- function call
|
| CFastFourierTransform::FFTR1DInternalEven(buf,m,buf3,plan);
|
| //--- function call
|
| CFastFourierTransform::FFTR1DInternalEven(buf2,m,buf3,plan);
|
| //--- change values
|
| buf[0]=buf[0]*buf2[0];
|
| buf[1]=buf[1]*buf2[1];
|
| //--- calculation
|
| for(i=1; i<=m/2-1; i++)
|
| {
|
| ax=buf[2*i+0];
|
| ay=buf[2*i+1];
|
| bx=buf2[2*i+0];
|
| by=buf2[2*i+1];
|
| tx=ax*bx-ay*by;
|
| ty=ax*by+ay*bx;
|
| buf[2*i+0]=tx;
|
| buf[2*i+1]=ty;
|
| }
|
| //--- function call
|
| CFastFourierTransform::FFTR1DInvInternalEven(buf,m,buf3,plan);
|
| //--- allocation
|
| ArrayResize(r,m);
|
| //--- copy
|
| for(i_=0; i_<=m-1; i_++)
|
| r[i_]=buf[i_];
|
| }
|
| else
|
| {
|
| //--- M is non-smooth or non-even,general code (circular/non-circular):
|
| //--- * first part is the same for circular and non-circular
|
| //--- convolutions. zero padding,FFTs,inverse FFTs
|
| //--- * second part differs:
|
| //--- * for non-circular convolution we just copy array
|
| //--- * for circular convolution we add array tail to its head
|
| p=CFtBase::FtBaseFindSmoothEven(m+n-1);
|
| //--- allocation
|
| ArrayResize(buf,p);
|
| for(i_=0; i_<=m-1; i_++)
|
| buf[i_]=a[i_];
|
| for(i=m; i<=p-1; i++)
|
| buf[i]=0;
|
| //--- allocation
|
| ArrayResize(buf2,p);
|
| for(i_=0; i_<=n-1; i_++)
|
| buf2[i_]=b[i_];
|
| for(i=n; i<=p-1; i++)
|
| buf2[i]=0;
|
| //--- allocation
|
| ArrayResize(buf3,p);
|
| //--- function call
|
| CFtBase::FtComplexFFTPlan(p/2,1,plan);
|
| //--- function call
|
| CFastFourierTransform::FFTR1DInternalEven(buf,p,buf3,plan);
|
| //--- function call
|
| CFastFourierTransform::FFTR1DInternalEven(buf2,p,buf3,plan);
|
| //--- change values
|
| buf[0]=buf[0]*buf2[0];
|
| buf[1]=buf[1]*buf2[1];
|
| //--- calculation
|
| for(i=1; i<=p/2-1; i++)
|
| {
|
| ax=buf[2*i+0];
|
| ay=buf[2*i+1];
|
| bx=buf2[2*i+0];
|
| by=buf2[2*i+1];
|
| tx=ax*bx-ay*by;
|
| ty=ax*by+ay*bx;
|
| buf[2*i+0]=tx;
|
| buf[2*i+1]=ty;
|
| }
|
| //--- function call
|
| CFastFourierTransform::FFTR1DInvInternalEven(buf,p,buf3,plan);
|
| //--- check
|
| if(circular)
|
| {
|
| //--- circular,add tail to head
|
| ArrayResize(r,m);
|
| for(i_=0; i_<=m-1; i_++)
|
| r[i_]=buf[i_];
|
| //--- check
|
| if(n>=2)
|
| {
|
| i1_=m;
|
| for(i_=0; i_<=n-2; i_++)
|
| r[i_]=r[i_]+buf[i_+i1_];
|
| }
|
| }
|
| else
|
| {
|
| //--- non-circular,just copy
|
| ArrayResize(r,m+n-1);
|
| for(i_=0; i_<=m+n-2; i_++)
|
| r[i_]=buf[i_];
|
| }
|
| }
|
| //--- exit the function
|
| return;
|
| }
|
| //--- overlap-add method
|
| if(alg==2)
|
| {
|
| //--- check
|
| if(!CAp::Assert((q+n-1)%2==0,__FUNCTION__+": internal error!"))
|
| return;
|
| //--- allocation
|
| ArrayResize(buf,q+n-1);
|
| ArrayResize(buf2,q+n-1);
|
| ArrayResize(buf3,q+n-1);
|
| //--- function call
|
| CFtBase::FtComplexFFTPlan((q+n-1)/2,1,plan);
|
| //--- prepare R
|
| if(circular)
|
| {
|
| //--- allocation
|
| ArrayResize(r,m);
|
| for(i=0; i<=m-1; i++)
|
| r[i]=0;
|
| }
|
| else
|
| {
|
| //--- allocation
|
| ArrayResize(r,m+n-1);
|
| for(i=0; i<=m+n-2; i++)
|
| r[i]=0;
|
| }
|
| //--- pre-calculated FFT(B)
|
| for(i_=0; i_<=n-1; i_++)
|
| buf2[i_]=b[i_];
|
| for(j=n; j<=q+n-2; j++)
|
| buf2[j]=0;
|
| //--- function call
|
| CFastFourierTransform::FFTR1DInternalEven(buf2,q+n-1,buf3,plan);
|
| //--- main overlap-add cycle
|
| i=0;
|
| //--- calculation
|
| while(i<=m-1)
|
| {
|
| p=MathMin(q,m-i);
|
| i1_=i;
|
| //--- copy
|
| for(i_=0; i_<=p-1; i_++)
|
| buf[i_]=a[i_+i1_];
|
| //--- make zero
|
| for(j=p; j<=q+n-2; j++)
|
| buf[j]=0;
|
| //--- function call
|
| CFastFourierTransform::FFTR1DInternalEven(buf,q+n-1,buf3,plan);
|
| //--- change values
|
| buf[0]=buf[0]*buf2[0];
|
| buf[1]=buf[1]*buf2[1];
|
| //--- calculation
|
| for(j=1; j<=(q+n-1)/2-1; j++)
|
| {
|
| ax=buf[2*j+0];
|
| ay=buf[2*j+1];
|
| bx=buf2[2*j+0];
|
| by=buf2[2*j+1];
|
| tx=ax*bx-ay*by;
|
| ty=ax*by+ay*bx;
|
| buf[2*j+0]=tx;
|
| buf[2*j+1]=ty;
|
| }
|
| //--- function call
|
| CFastFourierTransform::FFTR1DInvInternalEven(buf,q+n-1,buf3,plan);
|
| //--- check
|
| if(circular)
|
| {
|
| j1=MathMin(i+p+n-2,m-1)-i;
|
| j2=j1+1;
|
| }
|
| else
|
| {
|
| j1=p+n-2;
|
| j2=j1+1;
|
| }
|
| //--- change values
|
| i1_=-i;
|
| for(i_=i; i_<=i+j1; i_++)
|
| r[i_]=r[i_]+buf[i_+i1_];
|
| //--- check
|
| if(p+n-2>=j2)
|
| {
|
| i1_=j2;
|
| for(i_=0; i_<=p+n-2-j2; i_++)
|
| r[i_]=r[i_]+buf[i_+i1_];
|
| }
|
| i=i+p;
|
| }
|
| //--- exit the function
|
| return;
|
| }
|
| }
|
| //+------------------------------------------------------------------+
|
| //| Cross-correlation class |
|
| //+------------------------------------------------------------------+
|
| class CCorr
|
| {
|
| public:
|
| static void CorrC1D(complex &signal[],const int n,complex &pattern[],const int m,complex &r[]);
|
| static void CorrC1DCircular(complex &signal[],const int m,complex &pattern[],const int n,complex &c[]);
|
| static void CorrR1D(double &signal[],const int n,double &pattern[],const int m,double &r[]);
|
| static void CorrR1DCircular(double &signal[],const int m,double &pattern[],const int n,double &c[]);
|
| };
|
| //+------------------------------------------------------------------+
|
| //| 1-dimensional complex cross-correlation. |
|
| //| For given Pattern/Signal returns corr(Pattern,Signal) |
|
| //| (non-circular). |
|
| //| Correlation is calculated using reduction to convolution. |
|
| //| Algorithm with max(N,N)*log(max(N,N)) complexity is used (see |
|
| //| ConvC1D() for more info about performance). |
|
| //| IMPORTANT: |
|
| //| for historical reasons subroutine accepts its parameters in |
|
| //| reversed order: CorrC1D(Signal, Pattern) = Pattern x Signal |
|
| //| (using traditional definition of cross-correlation, denoting |
|
| //| cross-correlation as "x"). |
|
| //| INPUT PARAMETERS |
|
| //| Signal - array[0..N-1] - complex function to be |
|
| //| transformed, signal containing pattern |
|
| //| N - problem size |
|
| //| Pattern - array[0..M-1] - complex function to be |
|
| //| transformed, pattern to search withing signal |
|
| //| M - problem size |
|
| //| OUTPUT PARAMETERS |
|
| //| R - cross-correlation, array[0..N+M-2]: |
|
| //| * positive lags are stored in R[0..N-1], |
|
| //| R[i] = sum(conj(pattern[j])*signal[i+j] |
|
| //| * negative lags are stored in R[N..N+M-2], |
|
| //| R[N+M-1-i] = sum(conj(pattern[j])*signal[-i+j] |
|
| //| NOTE: |
|
| //| It is assumed that pattern domain is [0..M-1]. If Pattern is |
|
| //| non-zero on [-K..M-1], you can still use this subroutine, just |
|
| //| shift result by K. |
|
| //+------------------------------------------------------------------+
|
| void CCorr::CorrC1D(complex &signal[],const int n,complex &pattern[],
|
| const int m,complex &r[])
|
| {
|
| //--- create variables
|
| int i=0;
|
| int i_=0;
|
| int i1_=0;
|
| //--- create arrays
|
| complex p[];
|
| complex b[];
|
| //--- check
|
| if(!CAp::Assert(n>0 && m>0,__FUNCTION__+": incorrect N or M!"))
|
| return;
|
| //--- allocation
|
| ArrayResize(p,m);
|
| for(i=0; i<=m-1; i++)
|
| p[m-1-i]=CMath::Conj(pattern[i]);
|
| //--- function call
|
| CConv::ConvC1D(p,m,signal,n,b);
|
| //--- allocation
|
| ArrayResize(r,m+n-1);
|
| i1_=m-1;
|
| for(i_=0; i_<=n-1; i_++)
|
| r[i_]=b[i_+i1_];
|
| //--- check
|
| if(m+n-2>=n)
|
| {
|
| i1_=-n;
|
| for(i_=n; i_<=m+n-2; i_++)
|
| r[i_]=b[i_+i1_];
|
| }
|
| }
|
| //+------------------------------------------------------------------+
|
| //| 1-dimensional circular complex cross-correlation. |
|
| //| For given Pattern/Signal returns corr(Pattern,Signal) (circular).|
|
| //| Algorithm has linearithmic complexity for any M/N. |
|
| //| IMPORTANT: |
|
| //| for historical reasons subroutine accepts its parameters in |
|
| //| reversed order: CorrC1DCircular(Signal, Pattern) = Pattern x |
|
| //| Signal (using traditional definition of cross-correlation, |
|
| //| denoting cross-correlation as "x"). |
|
| //| INPUT PARAMETERS |
|
| //| Signal - array[0..N-1] - complex function to be |
|
| //| transformed, periodic signal containing pattern |
|
| //| N - problem size |
|
| //| Pattern - array[0..M-1] - complex function to be |
|
| //| transformed, non-periodic pattern to search |
|
| //| withing signal |
|
| //| M - problem size |
|
| //| OUTPUT PARAMETERS |
|
| //| R - convolution: A*B. array[0..M-1]. |
|
| //+------------------------------------------------------------------+
|
| void CCorr::CorrC1DCircular(complex &signal[],const int m,complex &pattern[],
|
| const int n,complex &c[])
|
| {
|
| //--- create variables
|
| int i1=0;
|
| int i2=0;
|
| int i=0;
|
| int j2=0;
|
| int i_=0;
|
| int i1_=0;
|
| //--- create arrays
|
| complex p[];
|
| complex b[];
|
| //--- check
|
| if(!CAp::Assert(n>0 && m>0,__FUNCTION__+": incorrect N or M!"))
|
| return;
|
| //--- normalize task: make M>=N,
|
| //--- so A will be longer (at least - not shorter) that B.
|
| if(m<n)
|
| {
|
| //--- allocation
|
| ArrayResize(b,m);
|
| //--- initialization
|
| for(i1=0; i1<=m-1; i1++)
|
| b[i1]=0;
|
| i1=0;
|
| //--- calculation
|
| while(i1<n)
|
| {
|
| //--- change values
|
| i2=MathMin(i1+m-1,n-1);
|
| j2=i2-i1;
|
| i1_=i1;
|
| for(i_=0; i_<=j2; i_++)
|
| b[i_]=b[i_]+pattern[i_+i1_];
|
| i1=i1+m;
|
| }
|
| //--- function call
|
| CorrC1DCircular(signal,m,b,m,c);
|
| //--- exit the function
|
| return;
|
| }
|
| //--- Task is normalized
|
| ArrayResize(p,n);
|
| for(i=0; i<=n-1; i++)
|
| p[n-1-i]=CMath::Conj(pattern[i]);
|
| //--- function call
|
| CConv::ConvC1DCircular(signal,m,p,n,b);
|
| //--- allocation
|
| ArrayResize(c,m);
|
| i1_=n-1;
|
| //--- copy
|
| for(i_=0; i_<=m-n; i_++)
|
| c[i_]=b[i_+i1_];
|
| //--- check
|
| if(m-n+1<=m-1)
|
| {
|
| i1_=-(m-n+1);
|
| for(i_=m-n+1; i_<=m-1; i_++)
|
| c[i_]=b[i_+i1_];
|
| }
|
| }
|
| //+------------------------------------------------------------------+
|
| //| 1-dimensional real cross-correlation. |
|
| //| For given Pattern/Signal returns corr(Pattern,Signal) |
|
| //| (non-circular). |
|
| //| Correlation is calculated using reduction to convolution. |
|
| //| Algorithm with max(N,N)*log(max(N,N)) complexity is used (see |
|
| //| ConvC1D() for more info about performance). |
|
| //| IMPORTANT: |
|
| //| for historical reasons subroutine accepts its parameters in |
|
| //| reversed order: CorrR1D(Signal, Pattern) = Pattern x Signal |
|
| //| (using traditional definition of cross-correlation, denoting|
|
| //| cross-correlation as "x"). |
|
| //| INPUT PARAMETERS |
|
| //| Signal - array[0..N-1] - real function to be transformed, |
|
| //| signal containing pattern |
|
| //| N - problem size |
|
| //| Pattern - array[0..M-1] - real function to be transformed, |
|
| //| pattern to search withing signal |
|
| //| M - problem size |
|
| //| OUTPUT PARAMETERS |
|
| //| R - cross-correlation, array[0..N+M-2]: |
|
| //| * positive lags are stored in R[0..N-1], |
|
| //| R[i] = sum(pattern[j]*signal[i+j] |
|
| //| * negative lags are stored in R[N..N+M-2], |
|
| //| R[N+M-1-i] = sum(pattern[j]*signal[-i+j] |
|
| //| NOTE: |
|
| //| It is assumed that pattern domain is [0..M-1]. If Pattern is |
|
| //| non-zero on [-K..M-1], you can still use this subroutine, just |
|
| //| shift result by K. |
|
| //+------------------------------------------------------------------+
|
| void CCorr::CorrR1D(double &signal[],const int n,double &pattern[],
|
| const int m,double &r[])
|
| {
|
| //--- create variables
|
| int i=0;
|
| int i_=0;
|
| int i1_=0;
|
| //--- create arrays
|
| double p[];
|
| double b[];
|
| //--- check
|
| if(!CAp::Assert(n>0 && m>0,__FUNCTION__+": incorrect N or M!"))
|
| return;
|
| //--- allocation
|
| ArrayResize(p,m);
|
| //--- copy
|
| for(i=0; i<=m-1; i++)
|
| p[m-1-i]=pattern[i];
|
| //--- function call
|
| CConv::ConvR1D(p,m,signal,n,b);
|
| //--- allocation
|
| ArrayResize(r,m+n-1);
|
| i1_=m-1;
|
| //--- copy
|
| for(i_=0; i_<=n-1; i_++)
|
| r[i_]=b[i_+i1_];
|
| //--- check
|
| if(m+n-2>=n)
|
| {
|
| i1_=-n;
|
| for(i_=n; i_<=m+n-2; i_++)
|
| r[i_]=b[i_+i1_];
|
| }
|
| }
|
| //+------------------------------------------------------------------+
|
| //| 1-dimensional circular real cross-correlation. |
|
| //| For given Pattern/Signal returns corr(Pattern,Signal) (circular).|
|
| //| Algorithm has linearithmic complexity for any M/N. |
|
| //| IMPORTANT: |
|
| //| for historical reasons subroutine accepts its parameters in |
|
| //| reversed order: CorrR1DCircular(Signal, Pattern) = Pattern x |
|
| //| Signal (using traditional definition of cross-correlation, |
|
| //| denoting cross-correlation as "x"). |
|
| //| INPUT PARAMETERS |
|
| //| Signal - array[0..N-1] - real function to be transformed, |
|
| //| periodic signal containing pattern |
|
| //| N - problem size |
|
| //| Pattern - array[0..M-1] - real function to be transformed, |
|
| //| non-periodic pattern to search withing signal |
|
| //| M - problem size |
|
| //| OUTPUT PARAMETERS |
|
| //| R - convolution: A*B. array[0..M-1]. |
|
| //+------------------------------------------------------------------+
|
| void CCorr::CorrR1DCircular(double &signal[],const int m,double &pattern[],
|
| const int n,double &c[])
|
| {
|
| //--- create variables
|
| int i1=0;
|
| int i2=0;
|
| int i=0;
|
| int j2=0;
|
| int i_=0;
|
| int i1_=0;
|
| //--- create arrays
|
| double p[];
|
| double b[];
|
| //--- check
|
| if(!CAp::Assert(n>0 && m>0,__FUNCTION__+": incorrect N or M!"))
|
| return;
|
| //--- normalize task: make M>=N,
|
| //--- so A will be longer (at least - not shorter) that B.
|
| if(m<n)
|
| {
|
| //--- allocation
|
| ArrayResize(b,m);
|
| //--- initialization
|
| for(i1=0; i1<=m-1; i1++)
|
| b[i1]=0;
|
| i1=0;
|
| //--- calculation
|
| while(i1<n)
|
| {
|
| //--- change values
|
| i2=MathMin(i1+m-1,n-1);
|
| j2=i2-i1;
|
| i1_=i1;
|
| for(i_=0; i_<=j2; i_++)
|
| b[i_]=b[i_]+pattern[i_+i1_];
|
| i1=i1+m;
|
| }
|
| //--- function call
|
| CorrR1DCircular(signal,m,b,m,c);
|
| //--- exit the function
|
| return;
|
| }
|
| //--- Task is normalized
|
| ArrayResize(p,n);
|
| for(i=0; i<=n-1; i++)
|
| p[n-1-i]=pattern[i];
|
| //--- function call
|
| CConv::ConvR1DCircular(signal,m,p,n,b);
|
| //--- allocation
|
| ArrayResize(c,m);
|
| i1_=n-1;
|
| //--- copy
|
| for(i_=0; i_<=m-n; i_++)
|
| c[i_]=b[i_+i1_];
|
| //--- check
|
| if(m-n+1<=m-1)
|
| {
|
| i1_=-(m-n+1);
|
| for(i_=m-n+1; i_<=m-1; i_++)
|
| c[i_]=b[i_+i1_];
|
| }
|
| }
|
| //+------------------------------------------------------------------+
|
| //| Fast Hartley Transform |
|
| //+------------------------------------------------------------------+
|
| class CFastHartleyTransform
|
| {
|
| public:
|
| static void FHTR1D(double &a[],const int n);
|
| static void FHTR1DInv(double &a[],const int n);
|
| };
|
| //+------------------------------------------------------------------+
|
| //| 1-dimensional Fast Hartley Transform. |
|
| //| Algorithm has O(N*logN) complexity for any N (composite or prime)|
|
| //| INPUT PARAMETERS |
|
| //| A - array[0..N-1] - real function to be transformed |
|
| //| N - problem size |
|
| //| OUTPUT PARAMETERS |
|
| //| A - FHT of a input array, array[0..N-1], |
|
| //| A_out[k]=sum(A_in[j]*(cos(2*pi*j*k/N)+sin(2*pi*j*k/N)),|
|
| //| j=0..N-1) |
|
| //+------------------------------------------------------------------+
|
| void CFastHartleyTransform::FHTR1D(double &a[],const int n)
|
| {
|
| //--- create a variable
|
| int i=0;
|
| //--- create array
|
| complex fa[];
|
| //--- object of class
|
| CFtPlan plan;
|
| //--- check
|
| if(!CAp::Assert(n>0,__FUNCTION__+": incorrect N!"))
|
| return;
|
| //--- Special case: N=1,FHT is just identity transform.
|
| //--- After this block we assume that N is strictly greater than 1.
|
| if(n==1)
|
| return;
|
| //--- Reduce FHt to real FFT
|
| CFastFourierTransform::FFTR1D(a,n,fa);
|
| //--- change values
|
| for(i=0; i<=n-1; i++)
|
| a[i]=fa[i].real-fa[i].imag;
|
| }
|
| //+------------------------------------------------------------------+
|
| //| 1-dimensional inverse FHT. |
|
| //| Algorithm has O(N*logN) complexity for any N (composite or prime)|
|
| //| INPUT PARAMETERS |
|
| //| A - array[0..N-1] - complex array to be transformed |
|
| //| N - problem size |
|
| //| OUTPUT PARAMETERS |
|
| //| A - inverse FHT of a input array, array[0..N-1] |
|
| //+------------------------------------------------------------------+
|
| void CFastHartleyTransform::FHTR1DInv(double &a[],const int n)
|
| {
|
| //--- create a variable
|
| int i=0;
|
| //--- check
|
| if(!CAp::Assert(n>0,__FUNCTION__+": incorrect N!"))
|
| return;
|
| //--- Special case: N=1,iFHT is just identity transform.
|
| //--- After this block we assume that N is strictly greater than 1.
|
| if(n==1)
|
| return;
|
| //--- Inverse FHT can be expressed in terms of the FHT as
|
| //--- invfht(x)=fht(x)/N
|
| FHTR1D(a,n);
|
| //--- change values
|
| for(i=0; i<=n-1; i++)
|
| a[i]=a[i]/n;
|
| }
|
| //+------------------------------------------------------------------+
|
