purebasic.info

PureBasic forum
Текущее время: Пн ноя 19, 2018 1:07 pm

Часовой пояс: UTC + 3 часа




Начать новую тему Ответить на тему  [ Сообщений: 5 ] 
Автор Сообщение
СообщениеДобавлено: Пн окт 15, 2018 7:03 pm 
Не в сети
профессор

Зарегистрирован: Чт сен 22, 2011 6:21 pm
Сообщений: 275
Благодарил (а): 36 раз.
Поблагодарили: 28 раз.
Пункты репутации: 0
Анализатор частоты на основе быстрого преобразования Фурье (модуль)
Обновление:
* исправлены длины при копировании массивов
зы. В тестовом блоке использован модуль линейного графика отсюда http://purebasic.info/phpBB3ex/viewtopic.php?f=10&t=4960
Код:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
;*********************************
; Module:   FFT
; Author:   void
; Compiler: PB546
; Ver:      002
;
;Анализатор частоты на основе быстрого преобразования Фурье
;Входной буффер (массив) должен быть размера кратного степени 2 (2,4,8,16,32,64,128...)
;Выходной буфер (массив) будет переопределен на размер половины входного
;Значения выходного буфера означают мощность текущей гармоники, а индекс означает количество периодов колебаний на кадр:
; индекс 0 - колебания с частотой 1 период на кадр или меньше
; индекс 1 - колебания с частотой 2 периоа на кадр
; индекс 2 - колебания с частотой 3 периоа на кадр
; верхняя граница частот ограничена частотой (размер входного кадра)/2 периодов.
; исходный алгоритм взят тут (вариант C++): "https://ru.wikibooks.org/wiki/Реализации_алгоритмов/Быстрое_преобразование_Фурье"
; оконные функции описаны тут: "https://ru.wikipedia.org/wiki/Оконное_преобразование_Фурье"
;*********************************
DeclareModule FFT
  EnableExplicit
 
  Declare Analyse(Array In.d(1), Array Ou.d(1))
  Declare Analyse2(Array In.d(2), Array Ou.d(2), Index=0)
  Declare WinFun(Array In.d(1), Array Ou.d(1), WinMode=1)
  Declare WinFun2(Array In.d(2), Array Ou.d(2), WinMode=1, Index=0)
 
EndDeclareModule
Module FFT
  EnableExplicit
  #Pi2 = #PI * 2.0
 
  Procedure WinFun(Array In.d(1), Array Ou.d(1), WinMode=1)
    Protected i, e = ArraySize(In()), N = e + 1
    Dim Ou(e)
    Select WinMode
      Case 0
        For i = 0 To e
          Ou(i) = In(i)
        Next
      Case 1  ;Окно Ханна
        For i = 0 To e
          Ou(i) = In(i) * (0.5 * (1 - Cos(#Pi2 * i / e)))
        Next
      Case 2  ;Окно Хемминга
        For i = 0 To e
          Ou(i) = In(i) * (0.53836 - 0.46164 * Cos(#Pi2 * i / e))
        Next
      Case 3  ;Окно Блэкмана
        For i = 0 To e
          Ou(i) = In(i) * (0.3635819 - 0.4891775 * Cos(#Pi2 * i / e) + 0.1365995 * Cos(2 * #Pi2 * i / e) - 0.0106411 * Cos(3 * #Pi2 * i / e))
        Next
    EndSelect
  EndProcedure
 
  Procedure WinFun2(Array In.d(2), Array Ou.d(2), WinMode=1, Index=0)
    Protected size = ArraySize(In(),2)
    Protected Dim In1.d(size), Dim Ou1.d(0)
    CopyMemory(@In(Index,0), @In1(0), (size + 1) * SizeOf(Double))
    WinFun(In1(), Ou1(), WinMode)
    Protected size1 = ArraySize(Ou(),1), size2 = ArraySize(Ou1())
    ReDim Ou(size1, size2)
    CopyMemory(@Ou1(), @Ou(Index,0), (size2 + 1) * SizeOf(Double))
    FreeArray(In1()): FreeArray(Ou1())
  EndProcedure
 
  Procedure Analyse(Array In.d(1), Array Ou.d(1))
    Protected.i i, j, n, m, Mmax, Istp, Nvl, Nft, MaxIndex
    Protected.d Tmpr, Tmpi, Wtmp, Theta
    Protected.d Wpr, Wpi, Wr, Wi, MaxValue
   
    Nvl = ArraySize(In()) + 1
    Nft = Nvl / 2
    Dim Ou.d(Nft - 1)
    n = Nvl * 2
    Dim Tmvl.d(n - 1)
   
    For i = 0 To n - 1 Step 2
      ;Tmvl(i) = 0 ;Тут нули после DIM
      Tmvl(i + 1) = In(i / 2)
    Next
   
    i = 1: j = 1
    While i < n
      If j > i
        Swap Tmvl(i), Tmvl(j)
        ;Swap Tmvl(i + 1), Tmvl(j + 1) ;Тут только нули
      EndIf
      i + 2: m = Nvl
      While (m >= 2) And (j > m)
        j - m: m >> 1
      Wend
      j + m
    Wend
   
    Mmax = 2
    While n > Mmax
      Theta = -#Pi2 / Mmax: Wpi = Sin(Theta)
      Wtmp = Sin(Theta / 2): Wpr = Wtmp * Wtmp * 2
      Istp = Mmax * 2: Wr = 1: Wi = 0: m = 1
     
      While m < Mmax
        i = m: m + 2: Tmpr = Wr: Tmpi = Wi
        Wr = Wr - Tmpr * Wpr - Tmpi * Wpi
        Wi = Wi + Tmpr * Wpi - Tmpi * Wpr
       
        While i < n
          j = i + Mmax
          Tmpr = Wr * Tmvl(j) - Wi * Tmvl(j - 1)
          Tmpi = Wi * Tmvl(j) + Wr * Tmvl(j - 1)
         
          Tmvl(j) = Tmvl(i) - Tmpr: Tmvl(j - 1) = Tmvl(i - 1) - Tmpi
          Tmvl(i) = Tmvl(i) + Tmpr: Tmvl(i - 1) = Tmvl(i - 1) + Tmpi
          i + Istp
        Wend
      Wend
     
      Mmax = Istp
    Wend
   
    For i = 0 To Nft - 1
      j = i * 2: Ou(i) = 2 * Sqr(Pow(Tmvl(j), 2) + Pow(Tmvl(j + 1), 2)) / Nvl
      If MaxValue < Ou(i)
        MaxValue = Ou(i)
        MaxIndex = i
      EndIf
    Next
   
    FreeArray(Tmvl())
    ProcedureReturn MaxIndex + 1
  EndProcedure
 
  Procedure Analyse2(Array In.d(2), Array Ou.d(2), Index=0)
    Protected size = ArraySize(In(),2), result
    Protected Dim In1.d(size), Dim Ou1.d(0)
    CopyMemory(@In(Index,0), @In1(0), (size + 1) * SizeOf(Double))
    result = Analyse(In1(), Ou1())
    Protected size1 = ArraySize(Ou(),1), size2 = ArraySize(Ou1())
    ReDim Ou(size1, size2)
    CopyMemory(@Ou1(), @Ou(Index,0), (size2 + 1) * SizeOf(Double))
    FreeArray(In1()): FreeArray(Ou1())
    ProcedureReturn result
  EndProcedure
 
EndModule
;Test Module
CompilerIf #PB_Compiler_IsMainFile
  EnableExplicit
 
  XIncludeFile "ChartMod.pbi"
 
  #WINSIDE = 1 << 10 ;должен быть кратен степени 2 (2,4,8,16,32,64...)
  Define i
 
  Dim d.d(0, #WINSIDE-1)
  Dim c.l(0, #WINSIDE-1)
  Dim d1.d(0,0)
  Dim d2.d(0,0)
 
 
  For i = 0 To ArraySize(d(), 2)
   
    d(0, i) = Sin(Radian(i * 10 * 360.0 / #WINSIDE)) +        ;сигнал частотой 10 периодов на кадр
              Sin(Radian(i * 250 * 360.0 / #WINSIDE)) * 1.1 + ;сигнал частотой 250 периодов на кадр
              Sin(Radian(i * 500 * 360.0 / #WINSIDE)) * 1.2    ;сигнал частотой 500 периодов на кадр
   
    c(0, i) = #Red
   
  Next
 
  CHART::Open("MultiSinus", 800, 620, d(), c(), #True)
 
  FFT::WinFun2(d(),d1(),1)
 
  CHART::Open("WinFunction", 800, 620, d1(), c(), #True)
 
  Debug "Максимальная гармоника (в периодах на кадр): " + Str(FFT::Analyse2(d1(), d2()))
 
  CHART::Open("FFT Analise", 800, 620, d2(), c(), #True)
 
CompilerEndIf
 



Последний раз редактировалось Kuzmat Вт окт 16, 2018 4:23 pm, всего редактировалось 2 раз(а).

Вернуться наверх
 Профиль  
 
СообщениеДобавлено: Пн окт 15, 2018 9:20 pm 
Не в сети
профессор

Зарегистрирован: Вс июл 05, 2009 5:55 pm
Сообщений: 310
Благодарил (а): 1 раз.
Поблагодарили: 10 раз.
Пункты репутации: 0
Kuzmat писал(а):
Анализатор частоты на основе быстрого преобразования Фурье (модуль)
зы. В тестовом блоке использован модуль линейного графика отсюда http://purebasic.info/phpBB3ex/viewtopic.php?f=10&t=4960

ругается на файлик :) которого нету, и тута и тама тоже: XIncludeFile "ChartMod.pbi".


Можно в двух словах про модули (я с таким неработал поэтому не понимаю), как они работают.
Код:
1
2
3
  CHART::Open("MultiSinus", 800, 620, d(), c(), #True)
  FFT::WinFun2(d(),d1(),1)
  CHART::Open("WinFunction", 800, 620, d1(), c(), #True)


_________________
искатель истины


Вернуться наверх
 Профиль  
 
СообщениеДобавлено: Пн окт 15, 2018 11:48 pm 
Не в сети
МОДЕРАТОР
Аватар пользователя

Зарегистрирован: Пн апр 09, 2007 4:53 pm
Сообщений: 11325
Благодарил (а): 4 раз.
Поблагодарили: 440 раз.
По ссылке есть код файла "ChartMod.pbi".

_________________
Компьютер позволяет решать все те проблемы, которые до его изобретения не существовали. :) :)


Вернуться наверх
 Профиль  
 
СообщениеДобавлено: Вт окт 16, 2018 8:16 am 
Не в сети
профессор

Зарегистрирован: Вс июл 05, 2009 5:55 pm
Сообщений: 310
Благодарил (а): 1 раз.
Поблагодарили: 10 раз.
Пункты репутации: 0
Пётр писал(а):
По ссылке есть код файла "ChartMod.pbi".

Ну ты и глазастый, ну у меня и браузер не находит такое-".pbi", где ты ее увидел?

_________________
искатель истины


Вернуться наверх
 Профиль  
 
СообщениеДобавлено: Вт окт 16, 2018 9:50 am 
Не в сети
профессор

Зарегистрирован: Пн июл 22, 2013 11:00 pm
Сообщений: 654
Благодарил (а): 2 раз.
Поблагодарили: 34 раз.
Пункты репутации: 9
balex1978 писал(а):
Пётр писал(а):
По ссылке есть код файла "ChartMod.pbi".

Ну ты и глазастый, ну у меня и браузер не находит такое-".pbi", где ты ее увидел?

Код сохрани под этим именем и расширением файла. :D


Вернуться наверх
 Профиль  
 
Показать сообщения за:  Сортировать по:  
Начать новую тему Ответить на тему  [ Сообщений: 5 ] 

Часовой пояс: UTC + 3 часа


Кто сейчас на форуме

Сейчас этот форум просматривают: нет зарегистрированных пользователей и гости: 4


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Перейти:  
Создано на основе phpBB® Forum Software © phpBB Group (блог о phpBB)
Сборка создана CMSart Studio
Русская поддержка phpBB