Visual Basic, .NET, ASP, VBScript
 

   
   
     

Форум - Работа с данными

Страница: 1 |

 

  Вопрос: алгоритм для БПФ Добавлено: 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
  1. ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
  2. 'Copyright (c) 2009, Sergey Bochkanov (ALGLIB project).
  3. '
  4. '>>> SOURCE LICENSE >>>
  5. 'This program is free software; you can redistribute it and/or modify
  6. 'it under the terms of the GNU General Public License as published by
  7. 'the Free Software Foundation (www.fsf.org); either version 2 of the
  8. 'License, or (at your option) any later version.
  9. '
  10. 'This program is distributed in the hope that it will be useful,
  11. 'but WITHOUT ANY WARRANTY; without even the implied warranty of
  12. 'MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13. 'GNU General Public License for more details.
  14. '
  15. 'A copy of the GNU General Public License is available at
  16. 'http://www.fsf.org/licensing/licenses
  17. '
  18. '>>> END OF LICENSE >>>
  19. ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
  20. 'Routines
  21. ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
  22. '1-dimensional complex FFT.
  23. '
  24. 'Array size N may be arbitrary number (composite or prime).  Composite  N's
  25. 'are handled with cache-oblivious variation of  a  Cooley-Tukey  algorithm.
  26. 'Small prime-factors are transformed using hard coded  codelets (similar to
  27. 'FFTW codelets, but without low-level  optimization),  large  prime-factors
  28. 'are handled with Bluestein's algorithm.
  29. '
  30. 'Fastests transforms are for smooth N's (prime factors are 2, 3,  5  only),
  31. 'most fast for powers of 2. When N have prime factors  larger  than  these,
  32. 'but orders of magnitude smaller than N, computations will be about 4 times
  33. 'slower than for nearby highly composite N's. When N itself is prime, speed
  34. 'will be 6 times lower.
  35. '
  36. 'Algorithm has O(N*logN) complexity for any N (composite or prime).
  37. '
  38. 'INPUT PARAMETERS
  39. '    A   -   array[0..N-1] - complex function to be transformed
  40. '    N   -   problem size
  41. '
  42. 'OUTPUT PARAMETERS
  43. '    A   -   DFT of a input array, array[0..N-1]
  44. '            A_out[j] = SUM(A_in[k]*exp(-2*pi*sqrt(-1)*j*k/N), k = 0..N-1)
  45. '
  46. '
  47. '  -- ALGLIB --
  48. '     Copyright 29.05.2009 by Bochkanov Sergey
  49. '
  50. ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
  51. Public Sub FFTC1D(ByRef A() As Complex, ByVal N As Long)
  52.     Dim Plan As FTPlan
  53.     Dim I As Long
  54.     Dim Buf() As Double
  55.  
  56.     
  57.     '
  58.     ' Special case: N=1, FFT is just identity transform.
  59.     ' After this block we assume that N is strictly greater than 1.
  60.     '
  61.     If N = 1# Then
  62.         Exit Sub
  63.     End If
  64.     
  65.     '
  66.     ' convert input array to the more convinient format
  67.     '
  68.     ReDim Buf(0 To 2# * N - 1)
  69.     For I = 0# To N - 1# Step 1
  70.         Buf(2# * I + 0#) = A(I).X
  71.         Buf(2# * I + 1#) = A(I).Y
  72.     Next I
  73.     
  74.     '
  75.     ' Generate plan and execute it.
  76.     '
  77.     ' Plan is a combination of a successive factorizations of N and
  78.     ' precomputed data. It is much like a FFTW plan, but is not stored
  79.     ' between subroutine calls and is much simpler.
  80.     '
  81.     Call FTBaseGenerateComplexFFTPlan(N, Plan)
  82.     Call FTBaseExecutePlan(Buf, 0#, N, Plan)
  83.     
  84.     '
  85.     ' result
  86.     '
  87.     For I = 0# To N - 1# Step 1
  88.         A(I).X = Buf(2# * I + 0#)
  89.         A(I).Y = Buf(2# * I + 1#)
  90.     Next I
  91. End Sub
  92.  
  93.  
  94. ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
  95. '1-dimensional complex inverse FFT.
  96. '
  97. 'Array size N may be arbitrary number (composite or prime).  Algorithm  has
  98. 'O(N*logN) complexity for any N (composite or prime).
  99. '
  100. 'See FFTC1D() description for more information about algorithm performance.
  101. '
  102. 'INPUT PARAMETERS
  103. '    A   -   array[0..N-1] - complex array to be transformed
  104. '    N   -   problem size
  105. '
  106. 'OUTPUT PARAMETERS
  107. '    A   -   inverse DFT of a input array, array[0..N-1]
  108. '            A_out[j] = SUM(A_in[k]/N*exp(+2*pi*sqrt(-1)*j*k/N), k = 0..N-1)
  109. '
  110. '
  111. '  -- ALGLIB --
  112. '     Copyright 29.05.2009 by Bochkanov Sergey
  113. '
  114. ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
  115. Public Sub FFTC1DInv(ByRef A() As Complex, ByVal N As Long)
  116.     Dim I As Long
  117.  
  118.     
  119.     '
  120.     ' Inverse DFT can be expressed in terms of the DFT as
  121.     '
  122.     '     invfft(x) = fft(x')'/N
  123.     '
  124.     ' here x' means conj(x).
  125.     '
  126.     For I = 0# To N - 1# Step 1
  127.         A(I).Y = -A(I).Y
  128.     Next I
  129.     Call FFTC1D(A, N)
  130.     For I = 0# To N - 1# Step 1
  131.         A(I).X = A(I).X / N
  132.         A(I).Y = -(A(I).Y / N)
  133.     Next I
  134. End Sub
  135.  
  136.  
  137. ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
  138. '1-dimensional real FFT.
  139. '
  140. 'Algorithm has O(N*logN) complexity for any N (composite or prime).
  141. '
  142. 'INPUT PARAMETERS
  143. '    A   -   array[0..N-1] - real function to be transformed
  144. '    N   -   problem size
  145. '
  146. 'OUTPUT PARAMETERS
  147. '    F   -   DFT of a input array, array[0..N-1]
  148. '            F[j] = SUM(A[k]*exp(-2*pi*sqrt(-1)*j*k/N), k = 0..N-1)
  149. '
  150. 'NOTE:
  151. '    F[] satisfies symmetry property F[k] = conj(F[N-k]),  so just one half
  152. 'of  array  is  usually needed. But for convinience subroutine returns full
  153. 'complex array (with frequencies above N/2), so its result may be  used  by
  154. 'other FFT-related subroutines.
  155. '
  156. '
  157. '  -- ALGLIB --
  158. '     Copyright 01.06.2009 by Bochkanov Sergey
  159. '
  160. ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
  161. Public Sub FFTR1D(ByRef A() As Double, ByVal N As Long, ByRef F() As Complex)
  162.     Dim I As Long
  163.     Dim N2 As Long
  164.     Dim Idx As Long
  165.     Dim Hn As Complex
  166.     Dim HmnC As Complex
  167.     Dim V As Complex
  168.     Dim Buf() As Double
  169.     Dim Plan As FTPlan
  170.     Dim i_ As Long
  171.  
  172.     
  173.     '
  174.     ' Special cases:
  175.     ' * N=1, FFT is just identity transform.
  176.     ' * N=2, FFT is simple too
  177.     '
  178.     ' After this block we assume that N is strictly greater than 2
  179.     '
  180.     If N = 1# Then
  181.         ReDim F(0 To 1# - 1)
  182.         F(0#) = C_Complex(A(0#))
  183.         Exit Sub
  184.     End If
  185.     If N = 2# Then
  186.         ReDim F(0 To 2# - 1)
  187.         F(0#).X = A(0#) + A(1#)
  188.         F(0#).Y = 0#
  189.         F(1#).X = A(0#) - A(1#)
  190.         F(1#).Y = 0#
  191.         Exit Sub
  192.     End If
  193.     
  194.     '
  195.     ' Choose between odd-size and even-size FFTs
  196.     '
  197.     If N Mod 2# = 0# Then
  198.         
  199.         '
  200.         ' even-size real FFT, use reduction to the complex task
  201.         '
  202.         N2 = N \ 2#
  203.         ReDim Buf(0 To N - 1)
  204.         For i_ = 0# To N - 1# Step 1
  205.             Buf(i_) = A(i_)
  206.         Next i_
  207.         Call FTBaseGenerateComplexFFTPlan(N2, Plan)
  208.         Call FTBaseExecutePlan(Buf, 0#, N2, Plan)
  209.         ReDim F(0 To N - 1)
  210.         For I = 0# To N2 Step 1
  211.             Idx = 2# * (I Mod N2)
  212.             Hn.X = Buf(Idx + 0#)
  213.             Hn.Y = Buf(Idx + 1#)
  214.             Idx = 2# * ((N2 - I) Mod N2)
  215.             HmnC.X = Buf(Idx + 0#)
  216.             HmnC.Y = -Buf(Idx + 1#)
  217.             V.X = -Sin(-(2# * Pi() * I / N))
  218.             V.Y = Cos(-(2# * Pi() * I / N))
  219.             F(I) = C_Sub(C_Add(Hn, HmnC), C_Mul(V, C_Sub(Hn, HmnC)))
  220.             F(I).X = 0.5 * F(I).X
  221.             F(I).Y = 0.5 * F(I).Y
  222.         Next I
  223.         For I = N2 + 1# To N - 1# Step 1
  224.             F(I) = Conj(F(N - I))
  225.         Next I
  226.         Exit Sub
  227.     Else
  228.         
  229.         '
  230.         ' use complex FFT
  231.         '
  232.         ReDim F(0 To N - 1)
  233.         For I = 0# To N - 1# Step 1
  234.             F(I) = C_Complex(A(I))
  235.         Next I
  236.         Call FFTC1D(F, N)
  237.         Exit Sub
  238.     End If
  239. End Sub
  240.  
  241.  
  242. ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
  243. '1-dimensional real inverse FFT.
  244. '
  245. 'Algorithm has O(N*logN) complexity for any N (composite or prime).
  246. '
  247. 'INPUT PARAMETERS
  248. '    F   -   array[0..floor(N/2)] - frequencies from forward real FFT
  249. '    N   -   problem size
  250. '
  251. 'OUTPUT PARAMETERS
  252. '    A   -   inverse DFT of a input array, array[0..N-1]
  253. '
  254. 'NOTE:
  255. '    F[] should satisfy symmetry property F[k] = conj(F[N-k]), so just  one
  256. 'half of frequencies array is needed - elements from 0 to floor(N/2).  F[0]
  257. 'is ALWAYS real. If N is even F[floor(N/2)] is real too. If N is odd,  then
  258. 'F[floor(N/2)] has no special properties.
  259. '
  260. 'Relying on properties noted above, FFTR1DInv subroutine uses only elements
  261. 'from 0th to floor(N/2)-th. It ignores imaginary part of F[0],  and in case
  262. 'N is even it ignores imaginary part of F[floor(N/2)] too.  So you can pass
  263. 'either frequencies array with N elements or reduced array with roughly N/2
  264. 'elements - subroutine will successfully transform both.
  265. '
  266. '
  267. '  -- ALGLIB --
  268. '     Copyright 01.06.2009 by Bochkanov Sergey
  269. '
  270. ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
  271. Public Sub FFTR1DInv(ByRef F() As Complex, _
  272.          ByVal N As Long, _
  273.          ByRef A() As Double)
  274.     Dim I As Long
  275.     Dim H() As Double
  276.     Dim FH() As Complex
  277.  
  278.     
  279.     '
  280.     ' Special case: N=1, FFT is just identity transform.
  281.     ' After this block we assume that N is strictly greater than 1.
  282.     '
  283.     If N = 1# Then
  284.         ReDim A(0 To 1# - 1)
  285.         A(0#) = F(0#).X
  286.         Exit Sub
  287.     End If
  288.     
  289.     '
  290.     ' inverse real FFT is reduced to the inverse real FHT,
  291.     ' which is reduced to the forward real FHT,
  292.     ' which is reduced to the forward real FFT.
  293.     '
  294.     ' Don't worry, it is really compact and efficient reduction :)
  295.     '
  296.     ReDim H(0 To N - 1)
  297.     ReDim A(0 To N - 1)
  298.     H(0#) = F(0#).X
  299.     For I = 1# To Int(N / 2#) - 1# Step 1
  300.         H(I) = F(I).X - F(I).Y
  301.         H(N - I) = F(I).X + F(I).Y
  302.     Next I
  303.     If N Mod 2# = 0# Then
  304.         H(Int(N / 2#)) = F(Int(N / 2#)).X
  305.     Else
  306.         H(Int(N / 2#)) = F(Int(N / 2#)).X - F(Int(N / 2#)).Y
  307.         H(Int(N / 2#) + 1#) = F(Int(N / 2#)).X + F(Int(N / 2#)).Y
  308.     End If
  309.     Call FFTR1D(H, N, FH)
  310.     For I = 0# To N - 1# Step 1
  311.         A(I) = (FH(I).X - FH(I).Y) / N
  312.     Next I
  313. End Sub
  314.  
  315.  
  316. ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
  317. 'Internal subroutine. Never call it directly!
  318. '
  319. '
  320. '  -- ALGLIB --
  321. '     Copyright 01.06.2009 by Bochkanov Sergey
  322. '
  323. ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
  324. Public Sub FFTR1DInternalEven(ByRef A() As Double, _
  325.          ByVal N As Long, _
  326.          ByRef Buf() As Double, _
  327.          ByRef Plan As FTPlan)
  328.     Dim X As Double
  329.     Dim Y As Double
  330.     Dim I As Long
  331.     Dim N2 As Long
  332.     Dim Idx As Long
  333.     Dim Hn As Complex
  334.     Dim HmnC As Complex
  335.     Dim V As Complex
  336.     Dim i_ As Long
  337.  
  338.     
  339.     '
  340.     ' Special cases:
  341.     ' * N=2
  342.     '
  343.     ' After this block we assume that N is strictly greater than 2
  344.     '
  345.     If N = 2# Then
  346.         X = A(0#) + A(1#)
  347.         Y = A(0#) - A(1#)
  348.         A(0#) = X
  349.         A(1#) = Y
  350.         Exit Sub
  351.     End If
  352.     
  353.     '
  354.     ' even-size real FFT, use reduction to the complex task
  355.     '
  356.     N2 = N \ 2#
  357.     For i_ = 0# To N - 1# Step 1
  358.         Buf(i_) = A(i_)
  359.     Next i_
  360.     Call FTBaseExecutePlan(Buf, 0#, N2, Plan)
  361.     A(0#) = Buf(0#) + Buf(1#)
  362.     For I = 1# To N2 - 1# Step 1
  363.         Idx = 2# * (I Mod N2)
  364.         Hn.X = Buf(Idx + 0#)
  365.         Hn.Y = Buf(Idx + 1#)
  366.         Idx = 2# * (N2 - I)
  367.         HmnC.X = Buf(Idx + 0#)
  368.         HmnC.Y = -Buf(Idx + 1#)
  369.         V.X = -Sin(-(2# * Pi() * I / N))
  370.         V.Y = Cos(-(2# * Pi() * I / N))
  371.         V = C_Sub(C_Add(Hn, HmnC), C_Mul(V, C_Sub(Hn, HmnC)))
  372.         A(2# * I + 0#) = 0.5 * V.X
  373.         A(2# * I + 1#) = 0.5 * V.Y
  374.     Next I
  375.     A(1#) = Buf(0#) - Buf(1#)
  376. End Sub
  377.  
  378.  
  379. ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
  380. 'Internal subroutine. Never call it directly!
  381. '
  382. '
  383. '  -- ALGLIB --
  384. '     Copyright 01.06.2009 by Bochkanov Sergey
  385. '
  386. ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
  387. Public Sub FFTR1DInvInternalEven(ByRef A() As Double, _
  388.          ByVal N As Long, _
  389.          ByRef Buf() As Double, _
  390.          ByRef Plan As FTPlan)
  391.     Dim X As Double
  392.     Dim Y As Double
  393.     Dim T As Double
  394.     Dim I As Long
  395.     Dim N2 As Long
  396.  
  397.     
  398.     '
  399.     ' Special cases:
  400.     ' * N=2
  401.     '
  402.     ' After this block we assume that N is strictly greater than 2
  403.     '
  404.     If N = 2# Then
  405.         X = 0.5 * (A(0#) + A(1#))
  406.         Y = 0.5 * (A(0#) - A(1#))
  407.         A(0#) = X
  408.         A(1#) = Y
  409.         Exit Sub
  410.     End If
  411.     
  412.     '
  413.     ' inverse real FFT is reduced to the inverse real FHT,
  414.     ' which is reduced to the forward real FHT,
  415.     ' which is reduced to the forward real FFT.
  416.     '
  417.     ' Don't worry, it is really compact and efficient reduction :)
  418.     '
  419.     N2 = N \ 2#
  420.     Buf(0#) = A(0#)
  421.     For I = 1# To N2 - 1# Step 1
  422.         X = A(2# * I + 0#)
  423.         Y = A(2# * I + 1#)
  424.         Buf(I) = X - Y
  425.         Buf(N - I) = X + Y
  426.     Next I
  427.     Buf(N2) = A(1#)
  428.     Call FFTR1DInternalEven(Buf, N, A, Plan)
  429.     A(0#) = Buf(0#) / N
  430.     T = 1# / N
  431.     For I = 1# To N2 - 1# Step 1
  432.         X = Buf(2# * I + 0#)
  433.         Y = Buf(2# * I + 1#)
  434.         A(I) = T * (X - Y)
  435.         A(N - I) = T * (X + Y)
  436.     Next I
  437.     A(N2) = Buf(1#) / N
  438. End Sub
  439.  
  440.  

Ответить

Номер ответа: 3
Автор ответа:
 AWP



ICQ: 345685652 

Вопросов: 96
Ответов: 1212
 Web-сайт: xawp.narod.ru
 Профиль | | #3
Добавлено: 01.03.13 01:05
Во! Это для FreeBasic писал, но смысл понятен должен быть.
  1.  
  2.     ' Áûñòðîå ïðåîáðàçîâàíèå Ôóðüå
  3.     Public Sub FastFourierTransform(ByVal a As Double PTR, ByVal nn As Long, ByVal InverseFFT As Integer) EXPORT
  4.         Dim ii As Long
  5.         Dim jj As Long
  6.         Dim n As Long
  7.         Dim mmax As Long
  8.         Dim m As Long
  9.         Dim j As Long
  10.         Dim istep As Long
  11.         Dim i As Long
  12.         Dim isign As Long
  13.         Dim wtemp As Double
  14.         Dim wr As Double
  15.         Dim wpr As Double
  16.         Dim wpi As Double
  17.         Dim wi As Double
  18.         Dim theta As Double
  19.         Dim tempr As Double
  20.         Dim tempi As Double
  21.  
  22.         If InverseFFT Then
  23.             isign = -1#
  24.         Else
  25.             isign = 1#
  26.         End If
  27.         n = 2# * nn
  28.         j = 1#
  29.         For ii = 1# To nn Step 1
  30.             i = 2# * ii - 1#
  31.             If j > i Then
  32.                 tempr = a[j - 1#]
  33.                 tempi = a[j]
  34.                 a[j - 1#] = a[i - 1#]
  35.                 a[j] = a
  36.                 a[i - 1#] = tempr
  37.                 a = tempi
  38.             End If
  39.             m = n \ 2#
  40.             Do While m >= 2# AND j > m
  41.                 j = j - m
  42.                 m = m \ 2#
  43.             Loop
  44.             j = j + m
  45.         Next ii
  46.         mmax = 2#
  47.         Do While n > mmax
  48.             istep = 2# * mmax
  49.             theta = 2# * 3.1415926535 / (isign * mmax)
  50.             wpr =  - (2# * dSQR(Sin(0.5 * theta)))
  51.             wpi = Sin(theta)
  52.             wr = 1#
  53.             wi = 0#
  54.             For ii = 1# To mmax \ 2# Step 1
  55.                 m = 2# * ii - 1#
  56.                 For jj = 0# To (n - m) \ istep Step 1
  57.                     i = m + jj * istep
  58.                     j = i + mmax
  59.                     tempr = wr * a[j - 1#] - wi * a[j]
  60.                     tempi = wr * a[j] + wi * a[j - 1#]
  61.                     a[j - 1#] = a[i - 1#] - tempr
  62.                     a[j] = a - tempi
  63.                     a[i - 1#] = a[i - 1#] + tempr
  64.                     a = a + tempi
  65.                 Next jj
  66.                 wtemp = wr
  67.                 wr = wr * wpr - wi * wpi + wr
  68.                 wi = wi * wpr + wtemp * wpi + wi
  69.             Next ii
  70.             mmax = istep
  71.         Loop
  72.         If InverseFFT Then
  73.             For i = 1# To 2# * nn Step 1
  74.                 a[i - 1#] = a[i - 1#] / nn
  75.             Next i
  76.         End If
  77.     End Sub

Ответить

Номер ответа: 4
Автор ответа:
 merr



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

Ответить

Страница: 1 |

Поиск по форуму





© Copyright 2002-2011 VBNet.RU | Пишите нам