VBNet - Форум - алгоритм для БПФ

Вопрос: алгоритм для БПФ Добавлено: 01.03.13 00:55

Автор вопроса: merr

Доброго времени суток, товарищи!
Может кто писал или видел пример алгоритма для быстрых преобразований Фурье (БПФ) на VB6?
Не имею морального права расспрашивать детали... Но буду признателен за ссылку ;)
Ответить

Ответы Всего ответов: 4

Номер ответа: 1 Автор ответа: AWP ICQ: 345685652 Вопросов: 96 Ответов: 1212	Web-сайт: xawp.narod.ru Профиль \| \| #1	Добавлено: 01.03.13 01:00
	Да есть. Как найду реализацию - скину. Ответить

Номер ответа: 2 Автор ответа: AWP ICQ: 345685652 Вопросов: 96 Ответов: 1212	Web-сайт: xawp.narod.ru Профиль \| \| #2	Добавлено: 01.03.13 01:02
	'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' 'Copyright (c) 2009, Sergey Bochkanov (ALGLIB project). ' '>>> SOURCE LICENSE >>> 'This program 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. ' 'A copy of the GNU General Public License is available at 'http://www.fsf.org/licensing/licenses ' '>>> END OF LICENSE >>> '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' 'Routines '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' '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(NlogN) 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(-2pisqrt(-1)jk/N), k = 0..N-1) ' ' ' -- ALGLIB -- ' Copyright 29.05.2009 by Bochkanov Sergey ' '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' Public Sub FFTC1D(ByRef A() As Complex, ByVal N As Long) Dim Plan As FTPlan Dim I As Long Dim Buf() As Double ' ' Special case: N=1, FFT is just identity transform. ' After this block we assume that N is strictly greater than 1. ' If N = 1# Then Exit Sub End If ' ' convert input array to the more convinient format ' ReDim Buf(0 To 2# * N - 1) For I = 0# To N - 1# Step 1 Buf(2# * I + 0#) = A(I).X Buf(2# * I + 1#) = A(I).Y Next I ' ' 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. ' Call FTBaseGenerateComplexFFTPlan(N, Plan) Call FTBaseExecutePlan(Buf, 0#, N, Plan) ' ' result ' For I = 0# To N - 1# Step 1 A(I).X = Buf(2# * I + 0#) A(I).Y = Buf(2# * I + 1#) Next I End Sub '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' '1-dimensional complex inverse FFT. ' 'Array size N may be arbitrary number (composite or prime). Algorithm has 'O(NlogN) 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]/Nexp(+2pisqrt(-1)jk/N), k = 0..N-1) ' ' ' -- ALGLIB -- ' Copyright 29.05.2009 by Bochkanov Sergey ' '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' Public Sub FFTC1DInv(ByRef A() As Complex, ByVal N As Long) Dim I As Long ' ' Inverse DFT can be expressed in terms of the DFT as ' ' invfft(x) = fft(x')'/N ' ' here x' means conj(x). ' For I = 0# To N - 1# Step 1 A(I).Y = -A(I).Y Next I Call FFTC1D(A, N) For I = 0# To N - 1# Step 1 A(I).X = A(I).X / N A(I).Y = -(A(I).Y / N) Next I End Sub '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' '1-dimensional real FFT. ' 'Algorithm has O(NlogN) 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(-2pisqrt(-1)jk/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. ' ' ' -- ALGLIB -- ' Copyright 01.06.2009 by Bochkanov Sergey ' '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' Public Sub FFTR1D(ByRef A() As Double, ByVal N As Long, ByRef F() As Complex) Dim I As Long Dim N2 As Long Dim Idx As Long Dim Hn As Complex Dim HmnC As Complex Dim V As Complex Dim Buf() As Double Dim Plan As FTPlan Dim i_ As Long ' ' 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# Then ReDim F(0 To 1# - 1) F(0#) = C_Complex(A(0#)) Exit Sub End If If N = 2# Then ReDim F(0 To 2# - 1) F(0#).X = A(0#) + A(1#) F(0#).Y = 0# F(1#).X = A(0#) - A(1#) F(1#).Y = 0# Exit Sub End If ' ' Choose between odd-size and even-size FFTs ' If N Mod 2# = 0# Then ' ' even-size real FFT, use reduction to the complex task ' N2 = N \ 2# ReDim Buf(0 To N - 1) For i_ = 0# To N - 1# Step 1 Buf(i_) = A(i_) Next i_ Call FTBaseGenerateComplexFFTPlan(N2, Plan) Call FTBaseExecutePlan(Buf, 0#, N2, Plan) ReDim F(0 To N - 1) For I = 0# To N2 Step 1 Idx = 2# * (I Mod N2) Hn.X = Buf(Idx + 0#) Hn.Y = Buf(Idx + 1#) Idx = 2# * ((N2 - I) Mod N2) HmnC.X = Buf(Idx + 0#) HmnC.Y = -Buf(Idx + 1#) V.X = -Sin(-(2# * Pi() * I / N)) V.Y = Cos(-(2# * Pi() * I / N)) F(I) = C_Sub(C_Add(Hn, HmnC), C_Mul(V, C_Sub(Hn, HmnC))) F(I).X = 0.5 * F(I).X F(I).Y = 0.5 * F(I).Y Next I For I = N2 + 1# To N - 1# Step 1 F(I) = Conj(F(N - I)) Next I Exit Sub Else ' ' use complex FFT ' ReDim F(0 To N - 1) For I = 0# To N - 1# Step 1 F(I) = C_Complex(A(I)) Next I Call FFTC1D(F, N) Exit Sub End If End Sub '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' '1-dimensional real inverse FFT. ' 'Algorithm has O(NlogN) 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. So you can pass 'either frequencies array with N elements or reduced array with roughly N/2 'elements - subroutine will successfully transform both. ' ' ' -- ALGLIB -- ' Copyright 01.06.2009 by Bochkanov Sergey ' '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' Public Sub FFTR1DInv(ByRef F() As Complex, _ ByVal N As Long, _ ByRef A() As Double) Dim I As Long Dim H() As Double Dim FH() As Complex ' ' Special case: N=1, FFT is just identity transform. ' After this block we assume that N is strictly greater than 1. ' If N = 1# Then ReDim A(0 To 1# - 1) A(0#) = F(0#).X Exit Sub End If ' ' 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 :) ' ReDim H(0 To N - 1) ReDim A(0 To N - 1) H(0#) = F(0#).X For I = 1# To Int(N / 2#) - 1# Step 1 H(I) = F(I).X - F(I).Y H(N - I) = F(I).X + F(I).Y Next I If N Mod 2# = 0# Then H(Int(N / 2#)) = F(Int(N / 2#)).X Else H(Int(N / 2#)) = F(Int(N / 2#)).X - F(Int(N / 2#)).Y H(Int(N / 2#) + 1#) = F(Int(N / 2#)).X + F(Int(N / 2#)).Y End If Call FFTR1D(H, N, FH) For I = 0# To N - 1# Step 1 A(I) = (FH(I).X - FH(I).Y) / N Next I End Sub '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' 'Internal subroutine. Never call it directly! ' ' ' -- ALGLIB -- ' Copyright 01.06.2009 by Bochkanov Sergey ' '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' Public Sub FFTR1DInternalEven(ByRef A() As Double, _ ByVal N As Long, _ ByRef Buf() As Double, _ ByRef Plan As FTPlan) Dim X As Double Dim Y As Double Dim I As Long Dim N2 As Long Dim Idx As Long Dim Hn As Complex Dim HmnC As Complex Dim V As Complex Dim i_ As Long ' ' Special cases: ' N=2 ' ' After this block we assume that N is strictly greater than 2 ' If N = 2# Then X = A(0#) + A(1#) Y = A(0#) - A(1#) A(0#) = X A(1#) = Y Exit Sub End If ' ' even-size real FFT, use reduction to the complex task ' N2 = N \ 2# For i_ = 0# To N - 1# Step 1 Buf(i_) = A(i_) Next i_ Call FTBaseExecutePlan(Buf, 0#, N2, Plan) A(0#) = Buf(0#) + Buf(1#) For I = 1# To N2 - 1# Step 1 Idx = 2# * (I Mod N2) Hn.X = Buf(Idx + 0#) Hn.Y = Buf(Idx + 1#) Idx = 2# * (N2 - I) HmnC.X = Buf(Idx + 0#) HmnC.Y = -Buf(Idx + 1#) V.X = -Sin(-(2# * Pi() * I / N)) V.Y = Cos(-(2# * Pi() * I / N)) V = C_Sub(C_Add(Hn, HmnC), C_Mul(V, C_Sub(Hn, HmnC))) A(2# * I + 0#) = 0.5 * V.X A(2# * I + 1#) = 0.5 * V.Y Next I A(1#) = Buf(0#) - Buf(1#) End Sub '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' 'Internal subroutine. Never call it directly! ' ' ' -- ALGLIB -- ' Copyright 01.06.2009 by Bochkanov Sergey ' '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' Public Sub FFTR1DInvInternalEven(ByRef A() As Double, _ ByVal N As Long, _ ByRef Buf() As Double, _ ByRef Plan As FTPlan) Dim X As Double Dim Y As Double Dim T As Double Dim I As Long Dim N2 As Long ' ' Special cases: ' * N=2 ' ' After this block we assume that N is strictly greater than 2 ' If N = 2# Then X = 0.5 * (A(0#) + A(1#)) Y = 0.5 * (A(0#) - A(1#)) A(0#) = X A(1#) = Y Exit Sub End If ' ' 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#) For I = 1# To N2 - 1# Step 1 X = A(2# * I + 0#) Y = A(2# * I + 1#) Buf(I) = X - Y Buf(N - I) = X + Y Next I Buf(N2) = A(1#) Call FFTR1DInternalEven(Buf, N, A, Plan) A(0#) = Buf(0#) / N T = 1# / N For I = 1# To N2 - 1# Step 1 X = Buf(2# * I + 0#) Y = Buf(2# * I + 1#) A(I) = T * (X - Y) A(N - I) = T * (X + Y) Next I A(N2) = Buf(1#) / N End Sub Ответить

Номер ответа: 3 Автор ответа: AWP ICQ: 345685652 Вопросов: 96 Ответов: 1212	Web-сайт: xawp.narod.ru Профиль \| \| #3	Добавлено: 01.03.13 01:05
	Во! Это для FreeBasic писал, но смысл понятен должен быть. ' Áûñòðîå ïðåîáðàçîâàíèå Ôóðüå Public Sub FastFourierTransform(ByVal a As Double PTR, ByVal nn As Long, ByVal InverseFFT As Integer) EXPORT Dim ii As Long Dim jj As Long Dim n As Long Dim mmax As Long Dim m As Long Dim j As Long Dim istep As Long Dim i As Long Dim isign As Long Dim wtemp As Double Dim wr As Double Dim wpr As Double Dim wpi As Double Dim wi As Double Dim theta As Double Dim tempr As Double Dim tempi As Double If InverseFFT Then isign = -1# Else isign = 1# End If n = 2# * nn j = 1# For ii = 1# To nn Step 1 i = 2# * ii - 1# If j > i Then tempr = a[j - 1#] tempi = a[j] a[j - 1#] = a[i - 1#] a[j] = a a[i - 1#] = tempr a = tempi End If m = n \ 2# Do While m >= 2# AND j > m j = j - m m = m \ 2# Loop j = j + m Next ii mmax = 2# Do While n > mmax istep = 2# * mmax theta = 2# * 3.1415926535 / (isign * mmax) wpr = - (2# * dSQR(Sin(0.5 * theta))) wpi = Sin(theta) wr = 1# wi = 0# For ii = 1# To mmax \ 2# Step 1 m = 2# * ii - 1# For jj = 0# To (n - m) \ istep Step 1 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 - tempi a[i - 1#] = a[i - 1#] + tempr a = a + tempi Next jj wtemp = wr wr = wr * wpr - wi * wpi + wr wi = wi * wpr + wtemp * wpi + wi Next ii mmax = istep Loop If InverseFFT Then For i = 1# To 2# * nn Step 1 a[i - 1#] = a[i - 1#] / nn Next i End If End Sub Ответить

Номер ответа: 4 Автор ответа: merr Вопросов: 11 Ответов: 31	Профиль \| \| #4	Добавлено: 01.03.13 01:22
	AWP, благодарю! Буду, так сказать, вникать в математику через Бейсик Ответить

Страница: 1 |