//+------------------------------------------------------------------+ //| dt_FFT.mq4 | //| Copyright © 2006, klot. | //| klot@mail.ru | //+------------------------------------------------------------------+ #property copyright "FFT for MQL4 klot" #property link "http://alglib.sources.ru/fft/" #property library #define pi 3.14159265358979323846 //+-------------------------------------------------------------------------------------------------------------------+ //| My function | //+-------------------------------------------------------------------------------------------------------------------+ //| #import "#_lib_FFT.ex4" | //| void fastfouriertransform(double& a[], int nn, bool inversefft); | //| void realfastfouriertransform(double& a[], int tnn, bool inversefft); | //| void fastcorellation(double& signal[], int signallen, double& pattern[], int patternlen); | //| void fastconvolution(double& signal[], int signallen, double& response[], int negativelen, int positivelen); | //| void fastsinetransform(double& a[], int tnn, bool inversefst); | //| void fastcosinetransform(double& a[], int tnn, bool inversefct); | //| void tworealffts(double a1[], double a2[], double& a[], double& b[], int tn); | //| #import | //+-------------------------------------------------------------------------------------------------------------------+ /************************************************************************* Быстрое преобразование Фурье Алгоритм проводит быстрое преобразование Фурье комплексной функции, заданной nn отсчетами на действительной оси. В зависимости от переданных параметров, может выполняться как прямое, так и обратное преобразование. Входные параметры: nn - Число значений функции. Должно быть степенью двойки. Алгоритм не проверяет правильность переданного значения. a - array [0 .. 2*nn-1] of Real Значения функции. I-ому значению соответствуют элементы a[2*I] (вещественная часть) и a[2*I+1] (мнимая часть). InverseFFT - направление преобразования. True, если обратное, False, если прямое. Выходные параметры: a - результат преобразования. Подробнее см. описание на сайте. http://alglib.sources.ru/fft/ *************************************************************************/ void fastfouriertransform(double& a[], int nn, bool inversefft) { int ii; int jj; int n; int mmax; int m; int j; int istep; int i; int isign; double wtemp; double wr; double wpr; double wpi; double wi; double theta; double tempr; double tempi; if( inversefft ) { isign = -1; } else { isign = 1; } n = 2*nn; j = 1; for(ii = 1; ii <= nn; ii++) { i = 2*ii-1; if( j>i ) { tempr = a[j-1]; tempi = a[j]; a[j-1] = a[i-1]; a[j] = a[i]; a[i-1] = tempr; a[i] = tempi; } m = n/2; while(m>=2&&j>m) { j = j-m; m = m/2; } j = j+m; } mmax = 2; while(n>mmax) { istep = 2*mmax; theta = 2.0*pi/(isign*mmax); wpr = -2.0*MathPow(MathSin(0.5*theta),2); wpi = MathSin(theta); wr = 1.0; wi = 0.0; for(ii = 1; ii <= mmax/2; ii++) { m = 2*ii-1; for(jj = 0; jj <= (n-m)/istep; jj++) { i = m+jj*istep; j = i+mmax; tempr = wr*a[j-1]-wi*a[j]; tempi = wr*a[j]+wi*a[j-1]; a[j-1] = a[i-1]-tempr; a[j] = a[i]-tempi; a[i-1] = a[i-1]+tempr; a[i] = a[i]+tempi; } wtemp = wr; wr = wr*wpr-wi*wpi+wr; wi = wi*wpr+wtemp*wpi+wi; } mmax = istep; } if( inversefft ) { for(i = 1; i <= 2*nn; i++) { a[i-1] = a[i-1]/nn; } } return; } /************************************************************************* Быстрое преобразование Фурье Алгоритм проводит быстрое преобразование Фурье вещественной функции, заданной n отсчетами на действительной оси. В зависимости от переданных параметров, может выполняться как прямое, так и обратное преобразование. Входные параметры: tnn - Число значений функции. Должно быть степенью двойки. Алгоритм не проверяет правильность переданного значения. a - array [0 .. nn-1] of Real Значения функции. InverseFFT - направление преобразования. True, если обратное, False, если прямое. Выходные параметры: a - результат преобразования. Подробнее см. описание на сайте. http://alglib.sources.ru/fft/ *************************************************************************/ void realfastfouriertransform(double& a[], int tnn, bool inversefft) { double twr; double twi; double twpr; double twpi; double twtemp; double ttheta; int i; int i1; int i2; int i3; int i4; double c1; double c2; double h1r; double h1i; double h2r; double h2i; double wrs; double wis; int nn; int ii; int jj; int n; int mmax; int m; int j; int istep; int isign; double wtemp; double wr; double wpr; double wpi; double wi; double theta; double tempr; double tempi; if( tnn==1 ) { return; } if( !inversefft ) { ttheta = 2.0*pi/tnn; c1 = 0.5; c2 = -0.5; } else { ttheta = 2.0*pi/tnn; c1 = 0.5; c2 = 0.5; ttheta = -ttheta; twpr = -2.0*MathPow(MathSin(0.5*ttheta),2); twpi = MathSin(ttheta); twr = 1.0+twpr; twi = twpi; for(i = 2; i <= tnn/4+1; i++) { i1 = i+i-2; i2 = i1+1; i3 = tnn+1-i2; i4 = i3+1; wrs = twr; wis = twi; h1r = c1*(a[i1]+a[i3]); h1i = c1*(a[i2]-a[i4]); h2r = -c2*(a[i2]+a[i4]); h2i = c2*(a[i1]-a[i3]); a[i1] = h1r+wrs*h2r-wis*h2i; a[i2] = h1i+wrs*h2i+wis*h2r; a[i3] = h1r-wrs*h2r+wis*h2i; a[i4] = -h1i+wrs*h2i+wis*h2r; twtemp = twr; twr = twr*twpr-twi*twpi+twr; twi = twi*twpr+twtemp*twpi+twi; } h1r = a[0]; a[0] = c1*(h1r+a[1]); a[1] = c1*(h1r-a[1]); } if( inversefft ) { isign = -1; } else { isign = 1; } n = tnn; nn = tnn/2; j = 1; for(ii = 1; ii <= nn; ii++) { i = 2*ii-1; if( j>i ) { tempr = a[j-1]; tempi = a[j]; a[j-1] = a[i-1]; a[j] = a[i]; a[i-1] = tempr; a[i] = tempi; } m = n/2; while(m>=2&&j>m) { j = j-m; m = m/2; } j = j+m; } mmax = 2; while(n>mmax) { istep = 2*mmax; theta = 2.0*pi/(isign*mmax); wpr = -2.0*MathPow(MathSin(0.5*theta),2); wpi = MathSin(theta); wr = 1.0; wi = 0.0; for(ii = 1; ii <= mmax/2; ii++) { m = 2*ii-1; for(jj = 0; jj <= (n-m)/istep; jj++) { i = m+jj*istep; j = i+mmax; tempr = wr*a[j-1]-wi*a[j]; tempi = wr*a[j]+wi*a[j-1]; a[j-1] = a[i-1]-tempr; a[j] = a[i]-tempi; a[i-1] = a[i-1]+tempr; a[i] = a[i]+tempi; } wtemp = wr; wr = wr*wpr-wi*wpi+wr; wi = wi*wpr+wtemp*wpi+wi; } mmax = istep; } if( inversefft ) { for(i = 1; i <= 2*nn; i++) { a[i-1] = a[i-1]/nn; } } if( !inversefft ) { twpr = -2.0*MathPow(MathSin(0.5*ttheta),2); twpi = MathSin(ttheta); twr = 1.0+twpr; twi = twpi; for(i = 2; i <= tnn/4+1; i++) { i1 = i+i-2; i2 = i1+1; i3 = tnn+1-i2; i4 = i3+1; wrs = twr; wis = twi; h1r = c1*(a[i1]+a[i3]); h1i = c1*(a[i2]-a[i4]); h2r = -c2*(a[i2]+a[i4]); h2i = c2*(a[i1]-a[i3]); a[i1] = h1r+wrs*h2r-wis*h2i; a[i2] = h1i+wrs*h2i+wis*h2r; a[i3] = h1r-wrs*h2r+wis*h2i; a[i4] = -h1i+wrs*h2i+wis*h2r; twtemp = twr; twr = twr*twpr-twi*twpi+twr; twi = twi*twpr+twtemp*twpi+twi; } h1r = a[0]; a[0] = h1r+a[1]; a[1] = h1r-a[1]; } return; } //=========================================================================================== /************************************************************************* Корелляция с использованием БПФ На входе: Signal - массив сигнал, с которым проводим корелляцию. Нумерация элементов от 0 до SignalLen-1 SignalLen - длина сигнала. Pattern - массив образец, корелляцию сигнала с которым мы ищем Нумерация элементов от 0 до PatternLen-1 PatternLen - длина образца На выходе: Signal - значения корелляции в точках от 0 до SignalLen-1. http://alglib.sources.ru/fft/ *************************************************************************/ void fastcorellation(double& signal[], int signallen, double& pattern[], int patternlen) { //ap::real_1d_array a1; //ap::real_1d_array a2; double a1[], a2[]; //--- int nl; int i; double t1; double t2; nl = signallen+patternlen; i = 1; while(ipositivelen ) { nl = nl+negativelen; } else { nl = nl+positivelen; } if( negativelen+1+positivelen>nl ) { nl = negativelen+1+positivelen; } i = 1; while(ii ) { tempr = a[j-1]; tempi = a[j]; a[j-1] = a[i-1]; a[j] = a[i]; a[i-1] = tempr; a[i] = tempi; } m = n/2; while(m>=2&&j>m) { j = j-m; m = m/2; } j = j+m; } mmax = 2; while(n>mmax) { istep = 2*mmax; theta = 2.0*pi/(isign*mmax); wpr = -2.0*MathPow(MathSin(0.5*theta),2); wpi = MathSin(theta); wr = 1.0; wi = 0.0; for(ii = 1; ii <= mmax/2; ii++) { m = 2*ii-1; for(jj = 0; jj <= (n-m)/istep; jj++) { i = m+jj*istep; j = i+mmax; tempr = wr*a[j-1]-wi*a[j]; tempi = wr*a[j]+wi*a[j-1]; a[j-1] = a[i-1]-tempr; a[j] = a[i]-tempi; a[i-1] = a[i-1]+tempr; a[i] = a[i]+tempi; } wtemp = wr; wr = wr*wpr-wi*wpi+wr; wi = wi*wpr+wtemp*wpi+wi; } mmax = istep; } twpr = -2.0*MathPow(MathSin(0.5*ttheta),2); twpi = MathSin(ttheta); twr = 1.0+twpr; twi = twpi; for(i = 2; i <= tnn/4+1; i++) { i1 = i+i-2; i2 = i1+1; i3 = tnn+1-i2; i4 = i3+1; wrs = twr; wis = twi; h1r = c1*(a[i1]+a[i3]); h1i = c1*(a[i2]-a[i4]); h2r = -c2*(a[i2]+a[i4]); h2i = c2*(a[i1]-a[i3]); a[i1] = h1r+wrs*h2r-wis*h2i; a[i2] = h1i+wrs*h2i+wis*h2r; a[i3] = h1r-wrs*h2r+wis*h2i; a[i4] = -h1i+wrs*h2i+wis*h2r; twtemp = twr; twr = twr*twpr-twi*twpi+twr; twi = twi*twpr+twtemp*twpi+twi; } h1r = a[0]; a[0] = h1r+a[1]; a[1] = h1r-a[1]; sum = 0.0; a[0] = 0.5*a[0]; a[1] = 0.0; for(jj = 0; jj <= tm-1; jj++) { j = 2*jj+1; sum = sum+a[j-1]; a[j-1] = a[j]; a[j] = sum; } if( inversefst ) { for(j = 1; j <= tnn; j++) { a[j-1] = a[j-1]*2/tnn; } } return; } //============================================================================================= /************************************************************************* Быстрое дискретное косинусное преобразование Алгоритм проводит быстрое косинусное преобразование вещественной функции, заданной nn отсчетами на действительной оси. В зависимости от переданных параметров, может выполняться как прямое, так и обратное преобразование. Входные параметры: tnn - Число значений функции минус один. Должно быть 1024 степенью двойки. Алгоритм не проверяет правильность переданного значения. a - array [0 .. nn] of Real 1025 Значения функции. InverseFCT - направление преобразования. True, если обратное, False, если прямое. Выходные параметры: a - результат преобразования. Подробнее см. описание на сайте. http://alglib.sources.ru/fft/ *************************************************************************/ void fastcosinetransform(double& a[], int tnn, bool inversefct) { int j; int n2; double sum; double y1; double y2; double theta; double wi; double wpi; double wr; double wpr; double wtemp; double twr; double twi; double twpr; double twpi; double twtemp; double ttheta; int i; int i1; int i2; int i3; int i4; double c1; double c2; double h1r; double h1i; double h2r; double h2i; double wrs; double wis; int nn; int ii; int jj; int n; int mmax; int m; int istep; int isign; double tempr; double tempi; if( tnn==1 ) { y1 = a[0]; y2 = a[1]; a[0] = 0.5*(y1+y2); a[1] = 0.5*(y1-y2); if( inversefct ) { a[0] = a[0]*2; a[1] = a[1]*2; } return; } wi = 0; wr = 1; theta = pi/tnn; wtemp = MathSin(theta*0.5); wpr = -2.0*wtemp*wtemp; wpi = MathSin(theta); sum = 0.5*(a[0]-a[tnn]); a[0] = 0.5*(a[0]+a[tnn]); n2 = tnn+2; for(j = 2; j <= tnn/2; j++) { wtemp = wr; wr = wtemp*wpr-wi*wpi+wtemp; wi = wi*wpr+wtemp*wpi+wi; y1 = 0.5*(a[j-1]+a[n2-j-1]); y2 = a[j-1]-a[n2-j-1]; a[j-1] = y1-wi*y2; a[n2-j-1] = y1+wi*y2; sum = sum+wr*y2; } ttheta = 2.0*pi/tnn; c1 = 0.5; c2 = -0.5; isign = 1; n = tnn; nn = tnn/2; j = 1; for(ii = 1; ii <= nn; ii++) { i = 2*ii-1; if( j>i ) { tempr = a[j-1]; tempi = a[j]; a[j-1] = a[i-1]; a[j] = a[i]; a[i-1] = tempr; a[i] = tempi; } m = n/2; while(m>=2&&j>m) { j = j-m; m = m/2; } j = j+m; } mmax = 2; while(n>mmax) { istep = 2*mmax; theta = 2.0*pi/(isign*mmax); wpr = -2.0*MathPow(MathSin(0.5*theta),2); wpi = MathSin(theta); wr = 1.0; wi = 0.0; for(ii = 1; ii <= mmax/2; ii++) { m = 2*ii-1; for(jj = 0; jj <= (n-m)/istep; jj++) { i = m+jj*istep; j = i+mmax; tempr = wr*a[j-1]-wi*a[j]; tempi = wr*a[j]+wi*a[j-1]; a[j-1] = a[i-1]-tempr; a[j] = a[i]-tempi; a[i-1] = a[i-1]+tempr; a[i] = a[i]+tempi; } wtemp = wr; wr = wr*wpr-wi*wpi+wr; wi = wi*wpr+wtemp*wpi+wi; } mmax = istep; } twpr = -2.0*MathPow(MathSin(0.5*ttheta),2); twpi = MathSin(ttheta); twr = 1.0+twpr; twi = twpi; for(i = 2; i <= tnn/4+1; i++) { i1 = i+i-2; i2 = i1+1; i3 = tnn+1-i2; i4 = i3+1; wrs = twr; wis = twi; h1r = c1*(a[i1]+a[i3]); h1i = c1*(a[i2]-a[i4]); h2r = -c2*(a[i2]+a[i4]); h2i = c2*(a[i1]-a[i3]); a[i1] = h1r+wrs*h2r-wis*h2i; a[i2] = h1i+wrs*h2i+wis*h2r; a[i3] = h1r-wrs*h2r+wis*h2i; a[i4] = -h1i+wrs*h2i+wis*h2r; twtemp = twr; twr = twr*twpr-twi*twpi+twr; twi = twi*twpr+twtemp*twpi+twi; } h1r = a[0]; a[0] = h1r+a[1]; a[1] = h1r-a[1]; a[tnn] = a[1]; a[1] = sum; j = 4; while(j<=tnn) { sum = sum+a[j-1]; a[j-1] = sum; j = j+2; } if( inversefct ) { for(j = 0; j <= tnn; j++) { a[j] = a[j]*2/tnn; } } return; } //=============================================================================================== /************************************************************************* Быстрое преобразование Фурье двух вещественных функций Алгоритм проводит быстрое преобразование Фурье двух вещественных функций, каждая из которых задана tn отсчетами на действительной оси. Алгоритм позволяет сэкономить время, но проводит только прямое преобразование. Входные параметры: tn - Число значений функций. Должно быть степенью двойки. Алгоритм не проверяет правильность переданного значения. a1 - array [0 .. nn-1] of Real Значения первой функции. a2 - array [0 .. nn-1] of Real Значения второй функции. Выходные параметры: a - Преобразование Фурье первой функции b - Преобразование Фурье второй функции (подробнее см. на сайте) http://alglib.sources.ru/fft/ *************************************************************************/ void tworealffts(double& a1[], double& a2[], double& a[], double& b[], int tn) { int jj; int j; double rep; double rem; double aip; double aim; int ii; int n; int nn; int mmax; int m; int istep; int i; int isign; double wtemp; double wr; double wpr; double wpi; double wi; double theta; double tempr; double tempi; nn = tn; //a.setbounds(0, 2*tn-1); //b.setbounds(0, 2*tn-1); ArrayResize(a, 2*tn); ArrayResize(b, 2*tn); //--- for(j = 1; j <= tn; j++) { jj = j+j; a[jj-2] = a1[j-1]; a[jj-1] = a2[j-1]; } isign = 1; n = 2*nn; j = 1; for(ii = 1; ii <= nn; ii++) { i = 2*ii-1; if( j>i ) { tempr = a[j-1]; tempi = a[j]; a[j-1] = a[i-1]; a[j] = a[i]; a[i-1] = tempr; a[i] = tempi; } m = n/2; while(m>=2&&j>m) { j = j-m; m = m/2; } j = j+m; } mmax = 2; while(n>mmax) { istep = 2*mmax; theta = 2.0*pi/(isign*mmax); wpr = -2.0*MathPow(MathSin(0.5*theta),2); wpi = MathSin(theta); wr = 1.0; wi = 0.0; for(ii = 1; ii <= mmax/2; ii++) { m = 2*ii-1; for(jj = 0; jj <= (n-m)/istep; jj++) { i = m+jj*istep; j = i+mmax; tempr = wr*a[j-1]-wi*a[j]; tempi = wr*a[j]+wi*a[j-1]; a[j-1] = a[i-1]-tempr; a[j] = a[i]-tempi; a[i-1] = a[i-1]+tempr; a[i] = a[i]+tempi; } wtemp = wr; wr = wr*wpr-wi*wpi+wr; wi = wi*wpr+wtemp*wpi+wi; } mmax = istep; } b[0] = a[1]; a[1] = 0.0; b[1] = 0.0; for(jj = 1; jj <= tn/2; jj++) { j = 2*jj+1; rep = 0.5*(a[j-1]+a[2*tn+1-j]); rem = 0.5*(a[j-1]-a[2*tn+1-j]); aip = 0.5*(a[j]+a[2*tn+2-j]); aim = 0.5*(a[j]-a[2*tn+2-j]); a[j-1] = rep; a[j] = aim; a[2*tn+1-j] = rep; a[2*tn+2-j] = -aim; b[j-1] = aip; b[j] = -rem; b[2*tn+1-j] = aip; b[2*tn+2-j] = rem; } return; } //==============================================================================================