|
| 1 | +package main |
| 2 | + |
| 3 | +import ( |
| 4 | + "fmt" |
| 5 | + "os" |
| 6 | + |
| 7 | + "github.com/go-audio/wav" |
| 8 | + "github.com/rawbits2010/LibBitDauer/package/assay" |
| 9 | +) |
| 10 | + |
| 11 | +func main() { |
| 12 | + |
| 13 | + file, err := os.Open("out.wav") |
| 14 | + if err != nil { |
| 15 | + panic(err) |
| 16 | + } |
| 17 | + defer file.Close() |
| 18 | + |
| 19 | + decoder := wav.NewDecoder(file) |
| 20 | + audioBuf, err := decoder.FullPCMBuffer() |
| 21 | + if err != nil { |
| 22 | + panic(err) |
| 23 | + } |
| 24 | + |
| 25 | + sampleRate := float64(audioBuf.Format.SampleRate) |
| 26 | + samples := audioBuf.AsFloatBuffer().Data |
| 27 | + |
| 28 | + windowedSignal := assay.ApplyHammingWindow(samples) |
| 29 | + spectrum := assay.GetFFTForFullSample(windowedSignal) |
| 30 | + |
| 31 | + frequencies, magnitudes, phases := assay.ExtractFFTResults(spectrum, sampleRate) |
| 32 | + |
| 33 | + dominantFreqIdx := assay.FindDominantFrequencyIdx(magnitudes) |
| 34 | + fmt.Printf("Dominant Frequency: %.2f Hz\n", frequencies[dominantFreqIdx]) |
| 35 | + |
| 36 | + err = assay.ExportFFTResultsToCSV(frequencies, magnitudes, phases, "fft_spectrum.csv") |
| 37 | + if err != nil { |
| 38 | + fmt.Println("Error writing file:", err) |
| 39 | + } else { |
| 40 | + fmt.Println("FFT Spectrum exported to fft_spectrum.csv") |
| 41 | + } |
| 42 | + |
| 43 | + err = assay.PlotMagnitudeAndPhase(frequencies, magnitudes, phases, "fft_spectrum.png") |
| 44 | + if err != nil { |
| 45 | + fmt.Println("Error writing file:", err) |
| 46 | + } else { |
| 47 | + fmt.Println("FFT Spectrum exported to fft_spectrum.png") |
| 48 | + } |
| 49 | + |
| 50 | + threshold := 10.0 |
| 51 | + peaks := assay.FindPeakIdxsWithPercentageThreshold(magnitudes, threshold) |
| 52 | + |
| 53 | + fmt.Println("Found Peaks:") |
| 54 | + for _, index := range peaks { |
| 55 | + fmt.Printf("Frequency: %.2f Hz, Magnitude: %.5f\n", frequencies[index], magnitudes[index]) |
| 56 | + } |
| 57 | + |
| 58 | + counts := assay.GetMagnitudePercentageDistribution(magnitudes, 100) |
| 59 | + err = assay.PlotBins(counts, 0, 100, 100, "fft_percdist.png") |
| 60 | + if err != nil { |
| 61 | + fmt.Println("Error writing file:", err) |
| 62 | + } else { |
| 63 | + fmt.Println("Magnitude distribution exported to fft_percdist.png") |
| 64 | + } |
| 65 | + |
| 66 | + waveform := IdentifyByHarmonicRatios(peaks, magnitudes) |
| 67 | + fmt.Printf("Identified waveform based on harmonics ratios: %s\n", waveform) |
| 68 | + |
| 69 | + waveform = IdentifyWaveformWithMetrics(peaks, magnitudes, frequencies) |
| 70 | + fmt.Printf("Identified waveform with metrics: %s\n", waveform) |
| 71 | +} |
| 72 | + |
| 73 | +// An example way of the identification process is to check the magnitude changes |
| 74 | +// in the harmonics. |
| 75 | +// This is only an example, peaks need to belong to a single signal. The ratios |
| 76 | +// and number of peaks are all tunable values, and they depend on the peak finding |
| 77 | +// method. |
| 78 | +func IdentifyByHarmonicRatios(peaks []int, magnitudes []float64) string { |
| 79 | + |
| 80 | + harmonicRatios := make([]float64, len(peaks)) |
| 81 | + for i, peakIdx := range peaks { |
| 82 | + harmonicRatios[i] = magnitudes[peakIdx] / magnitudes[0] |
| 83 | + } |
| 84 | + |
| 85 | + // values needs to be tuned |
| 86 | + if len(peaks) == 1 { |
| 87 | + // only base frequency |
| 88 | + return "Sine wave" |
| 89 | + } else if len(peaks) > 5 && (len(harmonicRatios) > 1 && harmonicRatios[1] > 0.5) { |
| 90 | + // base frequency + odd harmonics (3f, 5f, 7f...). |
| 91 | + return "Square wave" |
| 92 | + } else if len(peaks) <= 5 && (len(harmonicRatios) > 1 && harmonicRatios[1] < 0.5) { |
| 93 | + // base frequency + all odd harmonics, but amplitude shows heavily decay (e.g., 1/9, 1/25). |
| 94 | + return "Triangle wave" |
| 95 | + } else { |
| 96 | + // base frequency + all harmonics (odd and even) |
| 97 | + // or we have a complex waveform |
| 98 | + return "Sawtooth or Complex waveform" |
| 99 | + } |
| 100 | +} |
| 101 | + |
| 102 | +// An example way of the identification process is to check the bandwidth and |
| 103 | +// the spectral centroid of the signal. |
| 104 | +// This is only an example, peaks need to belong to a single signal. The ratios |
| 105 | +// and number of peaks are all tunable values, and they depend on the peak finding |
| 106 | +// method. |
| 107 | +func IdentifyWaveformWithMetrics(peaks []int, magnitudes, frequencies []float64) string { |
| 108 | + |
| 109 | + bandwidth := assay.CalculateBandwidth(peaks, frequencies) |
| 110 | + spectralCentroid := assay.CalculateSpectralCentroid(magnitudes, frequencies) |
| 111 | + |
| 112 | + // values needs to be tuned |
| 113 | + if bandwidth > 1000 && spectralCentroid > 1500 { |
| 114 | + // Wider bandwidth due to odd harmonics, with the spectral centroid shifted to higher frequencies. |
| 115 | + // High bandwidth and spectral centroid suggest a square wave. |
| 116 | + return "Square wave" |
| 117 | + } else if bandwidth < 300 && spectralCentroid < 1000 { |
| 118 | + // A sine wave has no harmonics, so its bandwidth is very small or zero. |
| 119 | + // The spectral centroid is at the fundamental frequency, as there are no higher harmonics. |
| 120 | + // If the bandwidth is small and the spectral centroid is at the fundamental, it's likely a sine wave. |
| 121 | + return "Sine wave" |
| 122 | + } else if bandwidth >= 300 && bandwidth <= 1000 { |
| 123 | + // Typically has a narrower bandwidth with harmonics decreasing with frequency. |
| 124 | + // The spectral centroid is lower than for square waves, as harmonics decrease relative to the fundamental. |
| 125 | + // If the bandwidth is medium and the spectral centroid is lower, it's likely a triangle wave. |
| 126 | + return "Triangle wave" |
| 127 | + } else { |
| 128 | + return "Complex waveform" |
| 129 | + } |
| 130 | +} |
0 commit comments