Вы не можете выбрать более 25 тем Темы должны начинаться с буквы или цифры, могут содержать дефисы(-) и должны содержать не более 35 символов.

430 строки
11KB

  1. package detector
  2. import (
  3. "math"
  4. "sort"
  5. "time"
  6. "sdr-visual-suite/internal/classifier"
  7. )
  8. type Event struct {
  9. ID int64 `json:"id"`
  10. Start time.Time `json:"start"`
  11. End time.Time `json:"end"`
  12. CenterHz float64 `json:"center_hz"`
  13. Bandwidth float64 `json:"bandwidth_hz"`
  14. PeakDb float64 `json:"peak_db"`
  15. SNRDb float64 `json:"snr_db"`
  16. FirstBin int `json:"first_bin"`
  17. LastBin int `json:"last_bin"`
  18. Class *classifier.Classification `json:"class,omitempty"`
  19. }
  20. type Detector struct {
  21. ThresholdDb float64
  22. MinDuration time.Duration
  23. Hold time.Duration
  24. EmaAlpha float64
  25. HysteresisDb float64
  26. MinStableFrames int
  27. GapTolerance time.Duration
  28. CFAREnabled bool
  29. CFARGuardCells int
  30. CFARTrainCells int
  31. CFARRank int
  32. CFARScaleDb float64
  33. binWidth float64
  34. nbins int
  35. sampleRate int
  36. ema []float64
  37. active map[int64]*activeEvent
  38. nextID int64
  39. }
  40. type activeEvent struct {
  41. id int64
  42. start time.Time
  43. lastSeen time.Time
  44. centerHz float64
  45. bwHz float64
  46. peakDb float64
  47. snrDb float64
  48. firstBin int
  49. lastBin int
  50. class *classifier.Classification
  51. stableHits int
  52. }
  53. type Signal struct {
  54. FirstBin int `json:"first_bin"`
  55. LastBin int `json:"last_bin"`
  56. CenterHz float64 `json:"center_hz"`
  57. BWHz float64 `json:"bw_hz"`
  58. PeakDb float64 `json:"peak_db"`
  59. SNRDb float64 `json:"snr_db"`
  60. Class *classifier.Classification `json:"class,omitempty"`
  61. }
  62. func New(thresholdDb float64, sampleRate int, fftSize int, minDur, hold time.Duration, emaAlpha, hysteresis float64, minStable int, gapTolerance time.Duration, cfarEnabled bool, cfarGuard, cfarTrain, cfarRank int, cfarScaleDb float64) *Detector {
  63. if minDur <= 0 {
  64. minDur = 250 * time.Millisecond
  65. }
  66. if hold <= 0 {
  67. hold = 500 * time.Millisecond
  68. }
  69. if emaAlpha <= 0 || emaAlpha > 1 {
  70. emaAlpha = 0.2
  71. }
  72. if hysteresis <= 0 {
  73. hysteresis = 3
  74. }
  75. if minStable <= 0 {
  76. minStable = 3
  77. }
  78. if gapTolerance <= 0 {
  79. gapTolerance = hold
  80. }
  81. if cfarGuard < 0 {
  82. cfarGuard = 2
  83. }
  84. if cfarTrain <= 0 {
  85. cfarTrain = 16
  86. }
  87. if cfarScaleDb <= 0 {
  88. cfarScaleDb = 6
  89. }
  90. if cfarRank <= 0 || cfarRank > 2*cfarTrain {
  91. cfarRank = int(math.Round(0.75 * float64(2*cfarTrain)))
  92. if cfarRank <= 0 {
  93. cfarRank = 1
  94. }
  95. }
  96. return &Detector{
  97. ThresholdDb: thresholdDb,
  98. MinDuration: minDur,
  99. Hold: hold,
  100. EmaAlpha: emaAlpha,
  101. HysteresisDb: hysteresis,
  102. MinStableFrames: minStable,
  103. GapTolerance: gapTolerance,
  104. CFAREnabled: cfarEnabled,
  105. CFARGuardCells: cfarGuard,
  106. CFARTrainCells: cfarTrain,
  107. CFARRank: cfarRank,
  108. CFARScaleDb: cfarScaleDb,
  109. binWidth: float64(sampleRate) / float64(fftSize),
  110. nbins: fftSize,
  111. sampleRate: sampleRate,
  112. ema: make([]float64, fftSize),
  113. active: map[int64]*activeEvent{},
  114. nextID: 1,
  115. }
  116. }
  117. func (d *Detector) Process(now time.Time, spectrum []float64, centerHz float64) ([]Event, []Signal) {
  118. signals := d.detectSignals(spectrum, centerHz)
  119. finished := d.matchSignals(now, signals)
  120. return finished, signals
  121. }
  122. // UpdateClasses refreshes active event classes from current signals.
  123. func (d *Detector) UpdateClasses(signals []Signal) {
  124. for _, s := range signals {
  125. for _, ev := range d.active {
  126. if overlapHz(s.CenterHz, s.BWHz, ev.centerHz, ev.bwHz) && math.Abs(s.CenterHz-ev.centerHz) < (s.BWHz+ev.bwHz)/2.0 {
  127. if s.Class != nil {
  128. if ev.class == nil || s.Class.Confidence >= ev.class.Confidence {
  129. ev.class = s.Class
  130. }
  131. }
  132. }
  133. }
  134. }
  135. }
  136. func (d *Detector) detectSignals(spectrum []float64, centerHz float64) []Signal {
  137. n := len(spectrum)
  138. if n == 0 {
  139. return nil
  140. }
  141. smooth := d.smoothSpectrum(spectrum)
  142. thresholds := d.cfarThresholds(smooth)
  143. noiseGlobal := median(smooth)
  144. var signals []Signal
  145. in := false
  146. start := 0
  147. peak := -1e9
  148. peakBin := 0
  149. for i := 0; i < n; i++ {
  150. v := smooth[i]
  151. thresholdOn := d.ThresholdDb
  152. if thresholds != nil && !math.IsNaN(thresholds[i]) {
  153. thresholdOn = thresholds[i]
  154. }
  155. thresholdOff := thresholdOn - d.HysteresisDb
  156. if v >= thresholdOn {
  157. if !in {
  158. in = true
  159. start = i
  160. peak = v
  161. peakBin = i
  162. } else if v > peak {
  163. peak = v
  164. peakBin = i
  165. }
  166. } else if in && v < thresholdOff {
  167. noise := noiseGlobal
  168. if thresholds != nil && peakBin >= 0 && peakBin < len(thresholds) && !math.IsNaN(thresholds[peakBin]) {
  169. noise = thresholds[peakBin] - d.CFARScaleDb
  170. }
  171. signals = append(signals, d.makeSignal(start, i-1, peak, peakBin, noise, centerHz, smooth))
  172. in = false
  173. }
  174. }
  175. if in {
  176. noise := noiseGlobal
  177. if thresholds != nil && peakBin >= 0 && peakBin < len(thresholds) && !math.IsNaN(thresholds[peakBin]) {
  178. noise = thresholds[peakBin] - d.CFARScaleDb
  179. }
  180. signals = append(signals, d.makeSignal(start, n-1, peak, peakBin, noise, centerHz, smooth))
  181. }
  182. return signals
  183. }
  184. func (d *Detector) makeSignal(first, last int, peak float64, peakBin int, noise float64, centerHz float64, spectrum []float64) Signal {
  185. centerBin := float64(first+last) / 2.0
  186. centerFreq := centerHz + (centerBin-float64(d.nbins)/2.0)*d.binWidth
  187. bw := float64(last-first+1) * d.binWidth
  188. snr := peak - noise
  189. return Signal{
  190. FirstBin: first,
  191. LastBin: last,
  192. CenterHz: centerFreq,
  193. BWHz: bw,
  194. PeakDb: peak,
  195. SNRDb: snr,
  196. }
  197. }
  198. func (d *Detector) cfarThresholds(spectrum []float64) []float64 {
  199. if !d.CFAREnabled || d.CFARTrainCells <= 0 {
  200. return nil
  201. }
  202. n := len(spectrum)
  203. train := d.CFARTrainCells
  204. guard := d.CFARGuardCells
  205. totalTrain := 2 * train
  206. if totalTrain <= 0 {
  207. return nil
  208. }
  209. rank := d.CFARRank
  210. if rank <= 0 || rank > totalTrain {
  211. rank = int(math.Round(0.75 * float64(totalTrain)))
  212. if rank <= 0 {
  213. rank = 1
  214. }
  215. }
  216. rankIdx := rank - 1
  217. thresholds := make([]float64, n)
  218. buf := make([]float64, 0, totalTrain)
  219. firstValid := guard + train
  220. lastValid := n - guard - train - 1
  221. for i := 0; i < n; i++ {
  222. if i < firstValid || i > lastValid {
  223. thresholds[i] = math.NaN()
  224. }
  225. }
  226. if firstValid > lastValid {
  227. return thresholds
  228. }
  229. // Build initial sorted window for first valid bin.
  230. leftStart := firstValid - guard - train
  231. buf = append(buf, spectrum[leftStart:leftStart+train]...)
  232. buf = append(buf, spectrum[firstValid+guard+1:firstValid+guard+1+train]...)
  233. sort.Float64s(buf)
  234. thresholds[firstValid] = buf[rankIdx] + d.CFARScaleDb
  235. // Slide window: remove outgoing bins and insert incoming bins (O(train) per step).
  236. for i := firstValid + 1; i <= lastValid; i++ {
  237. outLeft := spectrum[i-guard-train-1]
  238. outRight := spectrum[i+guard]
  239. inLeft := spectrum[i-guard-1]
  240. inRight := spectrum[i+guard+train]
  241. buf = removeValue(buf, outLeft)
  242. buf = removeValue(buf, outRight)
  243. buf = insertValue(buf, inLeft)
  244. buf = insertValue(buf, inRight)
  245. thresholds[i] = buf[rankIdx] + d.CFARScaleDb
  246. }
  247. return thresholds
  248. }
  249. func removeValue(sorted []float64, v float64) []float64 {
  250. if len(sorted) == 0 {
  251. return sorted
  252. }
  253. idx := sort.SearchFloat64s(sorted, v)
  254. if idx < len(sorted) && sorted[idx] == v {
  255. return append(sorted[:idx], sorted[idx+1:]...)
  256. }
  257. for i := idx - 1; i >= 0; i-- {
  258. if sorted[i] == v {
  259. return append(sorted[:i], sorted[i+1:]...)
  260. }
  261. if sorted[i] < v {
  262. break
  263. }
  264. }
  265. for i := idx + 1; i < len(sorted); i++ {
  266. if sorted[i] == v {
  267. return append(sorted[:i], sorted[i+1:]...)
  268. }
  269. if sorted[i] > v {
  270. break
  271. }
  272. }
  273. return sorted
  274. }
  275. func insertValue(sorted []float64, v float64) []float64 {
  276. idx := sort.SearchFloat64s(sorted, v)
  277. sorted = append(sorted, 0)
  278. copy(sorted[idx+1:], sorted[idx:])
  279. sorted[idx] = v
  280. return sorted
  281. }
  282. func (d *Detector) smoothSpectrum(spectrum []float64) []float64 {
  283. if d.ema == nil || len(d.ema) != len(spectrum) {
  284. d.ema = make([]float64, len(spectrum))
  285. copy(d.ema, spectrum)
  286. return d.ema
  287. }
  288. alpha := d.EmaAlpha
  289. for i := range spectrum {
  290. v := spectrum[i]
  291. d.ema[i] = alpha*v + (1-alpha)*d.ema[i]
  292. }
  293. // IMPORTANT: caller must not modify returned slice
  294. return d.ema
  295. }
  296. func (d *Detector) matchSignals(now time.Time, signals []Signal) []Event {
  297. used := make(map[int64]bool, len(d.active))
  298. for _, s := range signals {
  299. var best *activeEvent
  300. var candidates []struct {
  301. ev *activeEvent
  302. dist float64
  303. }
  304. for _, ev := range d.active {
  305. if overlapHz(s.CenterHz, s.BWHz, ev.centerHz, ev.bwHz) && math.Abs(s.CenterHz-ev.centerHz) < (s.BWHz+ev.bwHz)/2.0 {
  306. candidates = append(candidates, struct {
  307. ev *activeEvent
  308. dist float64
  309. }{ev: ev, dist: math.Abs(s.CenterHz - ev.centerHz)})
  310. }
  311. }
  312. if len(candidates) > 0 {
  313. sort.Slice(candidates, func(i, j int) bool { return candidates[i].dist < candidates[j].dist })
  314. best = candidates[0].ev
  315. }
  316. if best == nil {
  317. id := d.nextID
  318. d.nextID++
  319. d.active[id] = &activeEvent{
  320. id: id,
  321. start: now,
  322. lastSeen: now,
  323. centerHz: s.CenterHz,
  324. bwHz: s.BWHz,
  325. peakDb: s.PeakDb,
  326. snrDb: s.SNRDb,
  327. firstBin: s.FirstBin,
  328. lastBin: s.LastBin,
  329. class: s.Class,
  330. stableHits: 1,
  331. }
  332. continue
  333. }
  334. used[best.id] = true
  335. best.lastSeen = now
  336. best.stableHits++
  337. best.centerHz = (best.centerHz + s.CenterHz) / 2.0
  338. if s.BWHz > best.bwHz {
  339. best.bwHz = s.BWHz
  340. }
  341. if s.PeakDb > best.peakDb {
  342. best.peakDb = s.PeakDb
  343. }
  344. if s.SNRDb > best.snrDb {
  345. best.snrDb = s.SNRDb
  346. }
  347. if s.FirstBin < best.firstBin {
  348. best.firstBin = s.FirstBin
  349. }
  350. if s.LastBin > best.lastBin {
  351. best.lastBin = s.LastBin
  352. }
  353. if s.Class != nil {
  354. if best.class == nil || s.Class.Confidence >= best.class.Confidence {
  355. best.class = s.Class
  356. }
  357. }
  358. }
  359. var finished []Event
  360. for id, ev := range d.active {
  361. if used[id] {
  362. continue
  363. }
  364. if now.Sub(ev.lastSeen) < d.GapTolerance {
  365. continue
  366. }
  367. duration := ev.lastSeen.Sub(ev.start)
  368. if duration < d.MinDuration || ev.stableHits < d.MinStableFrames {
  369. delete(d.active, id)
  370. continue
  371. }
  372. finished = append(finished, Event{
  373. ID: ev.id,
  374. Start: ev.start,
  375. End: ev.lastSeen,
  376. CenterHz: ev.centerHz,
  377. Bandwidth: ev.bwHz,
  378. PeakDb: ev.peakDb,
  379. SNRDb: ev.snrDb,
  380. FirstBin: ev.firstBin,
  381. LastBin: ev.lastBin,
  382. Class: ev.class,
  383. })
  384. delete(d.active, id)
  385. }
  386. return finished
  387. }
  388. func overlapHz(c1, b1, c2, b2 float64) bool {
  389. l1 := c1 - b1/2.0
  390. r1 := c1 + b1/2.0
  391. l2 := c2 - b2/2.0
  392. r2 := c2 + b2/2.0
  393. return l1 <= r2 && l2 <= r1
  394. }
  395. func median(vals []float64) float64 {
  396. if len(vals) == 0 {
  397. return 0
  398. }
  399. cpy := append([]float64(nil), vals...)
  400. sort.Float64s(cpy)
  401. mid := len(cpy) / 2
  402. if len(cpy)%2 == 0 {
  403. return (cpy[mid-1] + cpy[mid]) / 2.0
  404. }
  405. return cpy[mid]
  406. }