You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

622 line
16KB

  1. package detector
  2. import (
  3. "math"
  4. "sort"
  5. "time"
  6. "sdr-visual-suite/internal/cfar"
  7. "sdr-visual-suite/internal/classifier"
  8. "sdr-visual-suite/internal/config"
  9. )
  10. type Event struct {
  11. ID int64 `json:"id"`
  12. Start time.Time `json:"start"`
  13. End time.Time `json:"end"`
  14. CenterHz float64 `json:"center_hz"`
  15. Bandwidth float64 `json:"bandwidth_hz"`
  16. PeakDb float64 `json:"peak_db"`
  17. SNRDb float64 `json:"snr_db"`
  18. FirstBin int `json:"first_bin"`
  19. LastBin int `json:"last_bin"`
  20. Class *classifier.Classification `json:"class,omitempty"`
  21. }
  22. type Detector struct {
  23. ThresholdDb float64
  24. MinDuration time.Duration
  25. Hold time.Duration
  26. EmaAlpha float64
  27. HysteresisDb float64
  28. MinStableFrames int
  29. GapTolerance time.Duration
  30. CFARScaleDb float64
  31. EdgeMarginDb float64
  32. MaxSignalBwHz float64
  33. MergeGapHz float64
  34. classHistorySize int
  35. classSwitchRatio float64
  36. binWidth float64
  37. nbins int
  38. sampleRate int
  39. ema []float64
  40. active map[int64]*activeEvent
  41. nextID int64
  42. cfarEngine cfar.CFAR
  43. lastThresholds []float64
  44. lastNoiseFloor float64
  45. }
  46. type activeEvent struct {
  47. id int64
  48. start time.Time
  49. lastSeen time.Time
  50. centerHz float64
  51. bwHz float64
  52. peakDb float64
  53. snrDb float64
  54. firstBin int
  55. lastBin int
  56. class *classifier.Classification
  57. stableHits int
  58. classHistory []classifier.SignalClass
  59. classIdx int
  60. }
  61. type Signal struct {
  62. FirstBin int `json:"first_bin"`
  63. LastBin int `json:"last_bin"`
  64. CenterHz float64 `json:"center_hz"`
  65. BWHz float64 `json:"bw_hz"`
  66. PeakDb float64 `json:"peak_db"`
  67. SNRDb float64 `json:"snr_db"`
  68. NoiseDb float64 `json:"noise_db,omitempty"`
  69. Class *classifier.Classification `json:"class,omitempty"`
  70. PLL *classifier.PLLResult `json:"pll,omitempty"`
  71. }
  72. func New(detCfg config.DetectorConfig, sampleRate int, fftSize int) *Detector {
  73. minDur := time.Duration(detCfg.MinDurationMs) * time.Millisecond
  74. hold := time.Duration(detCfg.HoldMs) * time.Millisecond
  75. gapTolerance := time.Duration(detCfg.GapToleranceMs) * time.Millisecond
  76. emaAlpha := detCfg.EmaAlpha
  77. hysteresis := detCfg.HysteresisDb
  78. minStable := detCfg.MinStableFrames
  79. cfarMode := detCfg.CFARMode
  80. binWidth := float64(sampleRate) / float64(fftSize)
  81. cfarGuard := int(math.Ceil(detCfg.CFARGuardHz / binWidth))
  82. if cfarGuard < 0 {
  83. cfarGuard = 0
  84. }
  85. cfarTrain := int(math.Ceil(detCfg.CFARTrainHz / binWidth))
  86. if cfarTrain < 1 {
  87. cfarTrain = 1
  88. }
  89. cfarRank := detCfg.CFARRank
  90. cfarScaleDb := detCfg.CFARScaleDb
  91. cfarWrap := detCfg.CFARWrapAround
  92. thresholdDb := detCfg.ThresholdDb
  93. edgeMarginDb := detCfg.EdgeMarginDb
  94. maxSignalBwHz := detCfg.MaxSignalBwHz
  95. mergeGapHz := detCfg.MergeGapHz
  96. classHistorySize := detCfg.ClassHistorySize
  97. classSwitchRatio := detCfg.ClassSwitchRatio
  98. if minDur <= 0 {
  99. minDur = 250 * time.Millisecond
  100. }
  101. if hold <= 0 {
  102. hold = 500 * time.Millisecond
  103. }
  104. if emaAlpha <= 0 || emaAlpha > 1 {
  105. emaAlpha = 0.2
  106. }
  107. if hysteresis <= 0 {
  108. hysteresis = 3
  109. }
  110. if minStable <= 0 {
  111. minStable = 3
  112. }
  113. if gapTolerance <= 0 {
  114. gapTolerance = hold
  115. }
  116. if cfarScaleDb <= 0 {
  117. cfarScaleDb = 6
  118. }
  119. if edgeMarginDb <= 0 {
  120. edgeMarginDb = 3.0
  121. }
  122. if maxSignalBwHz <= 0 {
  123. maxSignalBwHz = 150000
  124. }
  125. if mergeGapHz <= 0 {
  126. mergeGapHz = 5000
  127. }
  128. if classHistorySize <= 0 {
  129. classHistorySize = 10
  130. }
  131. if classSwitchRatio <= 0 || classSwitchRatio > 1 {
  132. classSwitchRatio = 0.6
  133. }
  134. if cfarRank <= 0 || cfarRank > 2*cfarTrain {
  135. cfarRank = int(math.Round(0.75 * float64(2*cfarTrain)))
  136. if cfarRank <= 0 {
  137. cfarRank = 1
  138. }
  139. }
  140. var cfarEngine cfar.CFAR
  141. if cfarMode != "" && cfarMode != "OFF" {
  142. cfarEngine = cfar.New(cfar.Config{
  143. Mode: cfar.Mode(cfarMode),
  144. GuardCells: cfarGuard,
  145. TrainCells: cfarTrain,
  146. Rank: cfarRank,
  147. ScaleDb: cfarScaleDb,
  148. WrapAround: cfarWrap,
  149. })
  150. }
  151. return &Detector{
  152. ThresholdDb: thresholdDb,
  153. MinDuration: minDur,
  154. Hold: hold,
  155. EmaAlpha: emaAlpha,
  156. HysteresisDb: hysteresis,
  157. MinStableFrames: minStable,
  158. GapTolerance: gapTolerance,
  159. CFARScaleDb: cfarScaleDb,
  160. EdgeMarginDb: edgeMarginDb,
  161. MaxSignalBwHz: maxSignalBwHz,
  162. MergeGapHz: mergeGapHz,
  163. classHistorySize: classHistorySize,
  164. classSwitchRatio: classSwitchRatio,
  165. binWidth: float64(sampleRate) / float64(fftSize),
  166. nbins: fftSize,
  167. sampleRate: sampleRate,
  168. ema: make([]float64, fftSize),
  169. active: map[int64]*activeEvent{},
  170. nextID: 1,
  171. cfarEngine: cfarEngine,
  172. }
  173. }
  174. func (d *Detector) Process(now time.Time, spectrum []float64, centerHz float64) ([]Event, []Signal) {
  175. signals := d.detectSignals(spectrum, centerHz)
  176. finished := d.matchSignals(now, signals)
  177. return finished, signals
  178. }
  179. func (d *Detector) LastThresholds() []float64 {
  180. if len(d.lastThresholds) == 0 {
  181. return nil
  182. }
  183. return append([]float64(nil), d.lastThresholds...)
  184. }
  185. func (d *Detector) LastNoiseFloor() float64 {
  186. return d.lastNoiseFloor
  187. }
  188. func (ev *activeEvent) updateClass(newCls *classifier.Classification, historySize int, switchRatio float64) {
  189. if newCls == nil {
  190. return
  191. }
  192. if historySize <= 0 {
  193. historySize = 10
  194. }
  195. if switchRatio <= 0 || switchRatio > 1 {
  196. switchRatio = 0.6
  197. }
  198. if len(ev.classHistory) != historySize {
  199. ev.classHistory = make([]classifier.SignalClass, historySize)
  200. ev.classIdx = 0
  201. }
  202. ev.classHistory[ev.classIdx%len(ev.classHistory)] = newCls.ModType
  203. ev.classIdx++
  204. if ev.class == nil {
  205. clone := *newCls
  206. ev.class = &clone
  207. return
  208. }
  209. counts := map[classifier.SignalClass]int{}
  210. filled := ev.classIdx
  211. if filled > len(ev.classHistory) {
  212. filled = len(ev.classHistory)
  213. }
  214. for i := 0; i < filled; i++ {
  215. c := ev.classHistory[i]
  216. if c != "" {
  217. counts[c]++
  218. }
  219. }
  220. var majority classifier.SignalClass
  221. majorityCount := 0
  222. for c, n := range counts {
  223. if n > majorityCount {
  224. majority = c
  225. majorityCount = n
  226. }
  227. }
  228. threshold := int(math.Ceil(float64(filled) * switchRatio))
  229. if threshold < 1 {
  230. threshold = 1
  231. }
  232. if majorityCount >= threshold && majority != ev.class.ModType {
  233. clone := *newCls
  234. clone.ModType = majority
  235. ev.class = &clone
  236. } else if majority == ev.class.ModType && newCls.Confidence > ev.class.Confidence {
  237. ev.class.Confidence = newCls.Confidence
  238. ev.class.Features = newCls.Features
  239. ev.class.SecondBest = newCls.SecondBest
  240. ev.class.Scores = newCls.Scores
  241. }
  242. }
  243. func (d *Detector) UpdateClasses(signals []Signal) {
  244. for _, s := range signals {
  245. for _, ev := range d.active {
  246. if overlapHz(s.CenterHz, s.BWHz, ev.centerHz, ev.bwHz) && math.Abs(s.CenterHz-ev.centerHz) < (s.BWHz+ev.bwHz)/2.0 {
  247. if s.Class != nil {
  248. ev.updateClass(s.Class, d.classHistorySize, d.classSwitchRatio)
  249. }
  250. }
  251. }
  252. }
  253. }
  254. func (d *Detector) detectSignals(spectrum []float64, centerHz float64) []Signal {
  255. n := len(spectrum)
  256. if n == 0 {
  257. return nil
  258. }
  259. smooth := d.smoothSpectrum(spectrum)
  260. var thresholds []float64
  261. if d.cfarEngine != nil {
  262. thresholds = d.cfarEngine.Thresholds(smooth)
  263. }
  264. d.lastThresholds = append(d.lastThresholds[:0], thresholds...)
  265. noiseGlobal := median(smooth)
  266. d.lastNoiseFloor = noiseGlobal
  267. var signals []Signal
  268. in := false
  269. start := 0
  270. peak := -1e9
  271. peakBin := 0
  272. for i := 0; i < n; i++ {
  273. v := smooth[i]
  274. thresholdOn := d.ThresholdDb
  275. if thresholds != nil && !math.IsNaN(thresholds[i]) {
  276. thresholdOn = thresholds[i]
  277. }
  278. thresholdOff := thresholdOn - d.HysteresisDb
  279. if v >= thresholdOn {
  280. if !in {
  281. in = true
  282. start = i
  283. peak = v
  284. peakBin = i
  285. } else if v > peak {
  286. peak = v
  287. peakBin = i
  288. }
  289. } else if in && v < thresholdOff {
  290. noise := noiseGlobal
  291. if thresholds != nil && peakBin >= 0 && peakBin < len(thresholds) && !math.IsNaN(thresholds[peakBin]) {
  292. noise = thresholds[peakBin] - d.CFARScaleDb
  293. }
  294. signals = append(signals, d.makeSignal(start, i-1, peak, peakBin, noise, centerHz, smooth))
  295. in = false
  296. }
  297. }
  298. if in {
  299. noise := noiseGlobal
  300. if thresholds != nil && peakBin >= 0 && peakBin < len(thresholds) && !math.IsNaN(thresholds[peakBin]) {
  301. noise = thresholds[peakBin] - d.CFARScaleDb
  302. }
  303. signals = append(signals, d.makeSignal(start, n-1, peak, peakBin, noise, centerHz, smooth))
  304. }
  305. signals = d.expandSignalEdges(signals, smooth, noiseGlobal, centerHz)
  306. for i := range signals {
  307. centerBin := float64(signals[i].FirstBin+signals[i].LastBin) / 2.0
  308. signals[i].CenterHz = d.centerFreqForBin(centerBin, centerHz)
  309. signals[i].BWHz = float64(signals[i].LastBin-signals[i].FirstBin+1) * d.binWidth
  310. }
  311. return signals
  312. }
  313. func (d *Detector) expandSignalEdges(signals []Signal, smooth []float64, noiseFloor float64, centerHz float64) []Signal {
  314. n := len(smooth)
  315. if n == 0 || len(signals) == 0 {
  316. return signals
  317. }
  318. margin := d.EdgeMarginDb
  319. if margin <= 0 {
  320. margin = 3.0
  321. }
  322. maxExpansionBins := int(d.MaxSignalBwHz / d.binWidth)
  323. if maxExpansionBins < 10 {
  324. maxExpansionBins = 10
  325. }
  326. for i := range signals {
  327. seed := signals[i]
  328. peakDb := seed.PeakDb
  329. localNoise := noiseFloor
  330. leftProbe := seed.FirstBin - 50
  331. rightProbe := seed.LastBin + 50
  332. if leftProbe >= 0 && rightProbe < n {
  333. leftNoise := minInRange(smooth, maxInt(0, leftProbe), maxInt(0, seed.FirstBin-5))
  334. rightNoise := minInRange(smooth, minInt(n-1, seed.LastBin+5), minInt(n-1, rightProbe))
  335. localNoise = math.Min(leftNoise, rightNoise)
  336. }
  337. edgeThreshold := localNoise + margin
  338. newFirst := seed.FirstBin
  339. prevVal := smooth[newFirst]
  340. for j := 0; j < maxExpansionBins; j++ {
  341. next := newFirst - 1
  342. if next < 0 {
  343. break
  344. }
  345. val := smooth[next]
  346. if val <= edgeThreshold {
  347. break
  348. }
  349. if val > prevVal+1.0 && val < peakDb-6.0 {
  350. break
  351. }
  352. prevVal = val
  353. newFirst = next
  354. }
  355. newLast := seed.LastBin
  356. prevVal = smooth[newLast]
  357. for j := 0; j < maxExpansionBins; j++ {
  358. next := newLast + 1
  359. if next >= n {
  360. break
  361. }
  362. val := smooth[next]
  363. if val <= edgeThreshold {
  364. break
  365. }
  366. if val > prevVal+1.0 && val < peakDb-6.0 {
  367. break
  368. }
  369. prevVal = val
  370. newLast = next
  371. }
  372. signals[i].FirstBin = newFirst
  373. signals[i].LastBin = newLast
  374. centerBin := float64(newFirst+newLast) / 2.0
  375. signals[i].CenterHz = d.centerFreqForBin(centerBin, centerHz)
  376. signals[i].BWHz = float64(newLast-newFirst+1) * d.binWidth
  377. }
  378. signals = d.mergeOverlapping(signals, centerHz)
  379. return signals
  380. }
  381. func (d *Detector) mergeOverlapping(signals []Signal, centerHz float64) []Signal {
  382. if len(signals) <= 1 {
  383. return signals
  384. }
  385. gapBins := 0
  386. if d.MergeGapHz > 0 && d.binWidth > 0 {
  387. gapBins = int(math.Ceil(d.MergeGapHz / d.binWidth))
  388. }
  389. sort.Slice(signals, func(i, j int) bool {
  390. return signals[i].FirstBin < signals[j].FirstBin
  391. })
  392. merged := []Signal{signals[0]}
  393. for i := 1; i < len(signals); i++ {
  394. last := &merged[len(merged)-1]
  395. cur := signals[i]
  396. gap := cur.FirstBin - last.LastBin - 1
  397. if gap <= gapBins {
  398. if cur.LastBin > last.LastBin {
  399. last.LastBin = cur.LastBin
  400. }
  401. if cur.PeakDb > last.PeakDb {
  402. last.PeakDb = cur.PeakDb
  403. }
  404. if cur.SNRDb > last.SNRDb {
  405. last.SNRDb = cur.SNRDb
  406. }
  407. centerBin := float64(last.FirstBin+last.LastBin) / 2.0
  408. last.BWHz = float64(last.LastBin-last.FirstBin+1) * d.binWidth
  409. last.CenterHz = d.centerFreqForBin(centerBin, centerHz)
  410. if cur.NoiseDb < last.NoiseDb || last.NoiseDb == 0 {
  411. last.NoiseDb = cur.NoiseDb
  412. }
  413. } else {
  414. merged = append(merged, cur)
  415. }
  416. }
  417. return merged
  418. }
  419. func (d *Detector) centerFreqForBin(bin float64, centerHz float64) float64 {
  420. return centerHz + (bin-float64(d.nbins)/2.0)*d.binWidth
  421. }
  422. func minInRange(s []float64, from, to int) float64 {
  423. if len(s) == 0 {
  424. return 0
  425. }
  426. if from < 0 {
  427. from = 0
  428. }
  429. if to >= len(s) {
  430. to = len(s) - 1
  431. }
  432. if from > to {
  433. return s[minInt(maxInt(from, 0), len(s)-1)]
  434. }
  435. m := s[from]
  436. for i := from + 1; i <= to; i++ {
  437. if s[i] < m {
  438. m = s[i]
  439. }
  440. }
  441. return m
  442. }
  443. func maxInt(a, b int) int {
  444. if a > b {
  445. return a
  446. }
  447. return b
  448. }
  449. func minInt(a, b int) int {
  450. if a < b {
  451. return a
  452. }
  453. return b
  454. }
  455. func overlapHz(center1, bw1, center2, bw2 float64) bool {
  456. left1 := center1 - bw1/2
  457. right1 := center1 + bw1/2
  458. left2 := center2 - bw2/2
  459. right2 := center2 + bw2/2
  460. return left1 <= right2 && left2 <= right1
  461. }
  462. func median(vals []float64) float64 {
  463. if len(vals) == 0 {
  464. return 0
  465. }
  466. cp := append([]float64(nil), vals...)
  467. sort.Float64s(cp)
  468. mid := len(cp) / 2
  469. if len(cp)%2 == 0 {
  470. return (cp[mid-1] + cp[mid]) / 2
  471. }
  472. return cp[mid]
  473. }
  474. func (d *Detector) makeSignal(first, last int, peak float64, peakBin int, noise float64, centerHz float64, spectrum []float64) Signal {
  475. centerBin := float64(first+last) / 2.0
  476. centerFreq := centerHz + (centerBin-float64(d.nbins)/2.0)*d.binWidth
  477. bw := float64(last-first+1) * d.binWidth
  478. snr := peak - noise
  479. return Signal{
  480. FirstBin: first,
  481. LastBin: last,
  482. CenterHz: centerFreq,
  483. BWHz: bw,
  484. PeakDb: peak,
  485. SNRDb: snr,
  486. NoiseDb: noise,
  487. }
  488. }
  489. func (d *Detector) smoothSpectrum(spectrum []float64) []float64 {
  490. if d.ema == nil || len(d.ema) != len(spectrum) {
  491. d.ema = make([]float64, len(spectrum))
  492. copy(d.ema, spectrum)
  493. return d.ema
  494. }
  495. alpha := d.EmaAlpha
  496. for i := range spectrum {
  497. v := spectrum[i]
  498. d.ema[i] = alpha*v + (1-alpha)*d.ema[i]
  499. }
  500. return d.ema
  501. }
  502. func (d *Detector) matchSignals(now time.Time, signals []Signal) []Event {
  503. used := make(map[int64]bool, len(d.active))
  504. for _, s := range signals {
  505. var best *activeEvent
  506. var candidates []struct {
  507. ev *activeEvent
  508. dist float64
  509. }
  510. for _, ev := range d.active {
  511. if overlapHz(s.CenterHz, s.BWHz, ev.centerHz, ev.bwHz) && math.Abs(s.CenterHz-ev.centerHz) < (s.BWHz+ev.bwHz)/2.0 {
  512. candidates = append(candidates, struct {
  513. ev *activeEvent
  514. dist float64
  515. }{ev: ev, dist: math.Abs(s.CenterHz - ev.centerHz)})
  516. }
  517. }
  518. if len(candidates) > 0 {
  519. sort.Slice(candidates, func(i, j int) bool { return candidates[i].dist < candidates[j].dist })
  520. best = candidates[0].ev
  521. }
  522. if best == nil {
  523. id := d.nextID
  524. d.nextID++
  525. d.active[id] = &activeEvent{
  526. id: id,
  527. start: now,
  528. lastSeen: now,
  529. centerHz: s.CenterHz,
  530. bwHz: s.BWHz,
  531. peakDb: s.PeakDb,
  532. snrDb: s.SNRDb,
  533. firstBin: s.FirstBin,
  534. lastBin: s.LastBin,
  535. class: s.Class,
  536. stableHits: 1,
  537. }
  538. continue
  539. }
  540. used[best.id] = true
  541. best.lastSeen = now
  542. best.stableHits++
  543. best.centerHz = (best.centerHz + s.CenterHz) / 2.0
  544. if best.bwHz <= 0 {
  545. best.bwHz = s.BWHz
  546. } else {
  547. const alpha = 0.15
  548. best.bwHz = alpha*s.BWHz + (1-alpha)*best.bwHz
  549. }
  550. if s.PeakDb > best.peakDb {
  551. best.peakDb = s.PeakDb
  552. }
  553. if s.SNRDb > best.snrDb {
  554. best.snrDb = s.SNRDb
  555. }
  556. if s.FirstBin < best.firstBin {
  557. best.firstBin = s.FirstBin
  558. }
  559. if s.LastBin > best.lastBin {
  560. best.lastBin = s.LastBin
  561. }
  562. if s.Class != nil {
  563. best.updateClass(s.Class, d.classHistorySize, d.classSwitchRatio)
  564. }
  565. }
  566. var finished []Event
  567. for id, ev := range d.active {
  568. if used[id] {
  569. continue
  570. }
  571. if now.Sub(ev.lastSeen) < d.GapTolerance {
  572. continue
  573. }
  574. duration := ev.lastSeen.Sub(ev.start)
  575. if duration < d.MinDuration || ev.stableHits < d.MinStableFrames {
  576. delete(d.active, id)
  577. continue
  578. }
  579. finished = append(finished, Event{
  580. ID: ev.id,
  581. Start: ev.start,
  582. End: ev.lastSeen,
  583. CenterHz: ev.centerHz,
  584. Bandwidth: ev.bwHz,
  585. PeakDb: ev.peakDb,
  586. SNRDb: ev.snrDb,
  587. FirstBin: ev.firstBin,
  588. LastBin: ev.lastBin,
  589. Class: ev.class,
  590. })
  591. delete(d.active, id)
  592. }
  593. return finished
  594. }