|
| 1 | +// Copyright 2026 Apple Inc. |
| 2 | +// |
| 3 | +// Use of this source code is governed by a BSD-3-clause license that can |
| 4 | +// be found in the LICENSE file or at https://opensource.org/licenses/BSD-3-Clause |
| 5 | + |
| 6 | +import Accelerate |
| 7 | +import AVFoundation |
| 8 | +import Foundation |
| 9 | + |
| 10 | +// Whisper mel spectrogram: sr=16000, n_fft=400, hop=160, n_mels=128 |
| 11 | +// Slaney-normalised filterbank, reflect-padded audio, matches WhisperFeatureExtractor. |
| 12 | +// |
| 13 | +// vDSP DFT only supports f×2^n sizes (f ∈ {1,3,5,15}); 400=5²×2⁴ doesn't qualify. |
| 14 | +// We precompute 201×400 DFT basis matrices and apply them with cblas_sgemv instead. |
| 15 | + |
| 16 | +enum WhisperMel { |
| 17 | + |
| 18 | + static let sampleRate: Double = 16_000 |
| 19 | + static let nFFT = 400 // analysis window (samples) |
| 20 | + static let hopLength = 160 |
| 21 | + static let nMelBins = 128 |
| 22 | + static let nFrames = 3_000 |
| 23 | + static let nSamples = 480_000 |
| 24 | + |
| 25 | + private static let nFreqs = nFFT / 2 + 1 // 201 |
| 26 | + |
| 27 | + // MARK: - Public |
| 28 | + |
| 29 | + static func fromFile(_ url: URL) throws -> [Float] { |
| 30 | + return fromPCM(try loadAndResample(url)) |
| 31 | + } |
| 32 | + |
| 33 | + // MARK: - Audio loading + resampling |
| 34 | + |
| 35 | + static func loadAndResample(_ url: URL) throws -> [Float] { |
| 36 | + let file = try AVAudioFile(forReading: url) |
| 37 | + let fmt = AVAudioFormat(commonFormat: .pcmFormatFloat32, |
| 38 | + sampleRate: sampleRate, channels: 1, interleaved: false)! |
| 39 | + guard let conv = AVAudioConverter(from: file.processingFormat, to: fmt) else { |
| 40 | + throw NSError(domain: "WhisperMel", code: 1, |
| 41 | + userInfo: [NSLocalizedDescriptionKey: |
| 42 | + "Cannot convert \(file.processingFormat) → 16 kHz mono"]) |
| 43 | + } |
| 44 | + let cap = AVAudioFrameCount( |
| 45 | + ceil(Double(file.length) * sampleRate / file.processingFormat.sampleRate) + 1) |
| 46 | + let out = AVAudioPCMBuffer(pcmFormat: fmt, frameCapacity: cap)! |
| 47 | + var fed = false; var err: NSError? |
| 48 | + conv.convert(to: out, error: &err) { _, status in |
| 49 | + guard !fed else { status.pointee = .endOfStream; return nil } |
| 50 | + fed = true |
| 51 | + let buf = AVAudioPCMBuffer(pcmFormat: file.processingFormat, |
| 52 | + frameCapacity: AVAudioFrameCount(file.length))! |
| 53 | + try? file.read(into: buf) |
| 54 | + status.pointee = buf.frameLength > 0 ? .haveData : .endOfStream |
| 55 | + return buf |
| 56 | + } |
| 57 | + if let e = err { throw e } |
| 58 | + return Array(UnsafeBufferPointer(start: out.floatChannelData![0], |
| 59 | + count: Int(out.frameLength))) |
| 60 | + } |
| 61 | + |
| 62 | + // MARK: - Precomputed DFT basis (201 × 400) |
| 63 | + // cos_basis[k, n] = cos(2π k n / 400) → Y[k].real = cos_basis @ x |
| 64 | + // sin_basis[k, n] = -sin(2π k n / 400) → Y[k].imag = sin_basis @ x |
| 65 | + |
| 66 | + static let cosBasis: [Float] = { |
| 67 | + var m = [Float](repeating: 0, count: (nFFT / 2 + 1) * nFFT) |
| 68 | + for k in 0...nFFT / 2 { |
| 69 | + for n in 0..<nFFT { |
| 70 | + m[k * nFFT + n] = cos(2 * Float.pi * Float(k) * Float(n) / Float(nFFT)) |
| 71 | + } |
| 72 | + } |
| 73 | + return m |
| 74 | + }() |
| 75 | + |
| 76 | + static let sinBasis: [Float] = { |
| 77 | + var m = [Float](repeating: 0, count: (nFFT / 2 + 1) * nFFT) |
| 78 | + for k in 0...nFFT / 2 { |
| 79 | + for n in 0..<nFFT { |
| 80 | + m[k * nFFT + n] = -sin(2 * Float.pi * Float(k) * Float(n) / Float(nFFT)) |
| 81 | + } |
| 82 | + } |
| 83 | + return m |
| 84 | + }() |
| 85 | + |
| 86 | + // MARK: - Mel filterbank (128 × 201, Slaney-normalised) |
| 87 | + |
| 88 | + static let melFilterbank: [Float] = makeMelFilterbank() |
| 89 | + |
| 90 | + // MARK: - Mel computation |
| 91 | + |
| 92 | + static func fromPCM(_ raw: [Float]) -> [Float] { |
| 93 | + // 1. Trim / zero-pad to nSamples |
| 94 | + var audio = raw |
| 95 | + if audio.count > nSamples { audio = Array(audio.prefix(nSamples)) } |
| 96 | + else if audio.count < nSamples { |
| 97 | + audio += [Float](repeating: 0, count: nSamples - audio.count) |
| 98 | + } |
| 99 | + |
| 100 | + // 2. Reflect-pad by nFFT/2 (matches np.pad(..., mode='reflect')) |
| 101 | + let pad = nFFT / 2 // 200 |
| 102 | + var padded = [Float](repeating: 0, count: nSamples + 2 * pad) |
| 103 | + for i in 0..<pad { padded[pad - 1 - i] = audio[i + 1] } |
| 104 | + for i in 0..<nSamples { padded[pad + i] = audio[i] } |
| 105 | + for i in 0..<pad { padded[pad + nSamples + i] = audio[nSamples - 2 - i] } |
| 106 | + |
| 107 | + // 3. Hann window |
| 108 | + var window = [Float](repeating: 0, count: nFFT) |
| 109 | + for i in 0..<nFFT { |
| 110 | + window[i] = Float(0.5 * (1 - cos(2 * Double.pi * Double(i) / Double(nFFT - 1)))) |
| 111 | + } |
| 112 | + |
| 113 | + var frame = [Float](repeating: 0, count: nFFT) |
| 114 | + var yReal = [Float](repeating: 0, count: nFreqs) |
| 115 | + var yImag = [Float](repeating: 0, count: nFreqs) |
| 116 | + var powerSpec = [Float](repeating: 0, count: nFreqs) |
| 117 | + var melFrame = [Float](repeating: 0, count: nMelBins) |
| 118 | + var mel = [Float](repeating: 0, count: nMelBins * nFrames) |
| 119 | + |
| 120 | + for t in 0..<nFrames { |
| 121 | + let offset = t * hopLength |
| 122 | + |
| 123 | + // Apply Hann window |
| 124 | + vDSP_vmul(Array(padded[offset ..< offset + nFFT]), 1, |
| 125 | + window, 1, &frame, 1, vDSP_Length(nFFT)) |
| 126 | + |
| 127 | + // DFT via matrix multiply: Y[k] = cosBasis[k,:] @ frame - i × sinBasis[k,:] @ frame |
| 128 | + cblas_sgemv(CblasRowMajor, CblasNoTrans, |
| 129 | + Int32(nFreqs), Int32(nFFT), 1.0, cosBasis, Int32(nFFT), |
| 130 | + frame, 1, 0.0, &yReal, 1) |
| 131 | + cblas_sgemv(CblasRowMajor, CblasNoTrans, |
| 132 | + Int32(nFreqs), Int32(nFFT), 1.0, sinBasis, Int32(nFFT), |
| 133 | + frame, 1, 0.0, &yImag, 1) |
| 134 | + |
| 135 | + // Power spectrum |Y[k]|² = yReal² + yImag² |
| 136 | + vDSP_vmma(yReal, 1, yReal, 1, yImag, 1, yImag, 1, &powerSpec, 1, vDSP_Length(nFreqs)) |
| 137 | + |
| 138 | + // Apply mel filterbank: (128×201) × (201) → (128) |
| 139 | + cblas_sgemv(CblasRowMajor, CblasNoTrans, |
| 140 | + Int32(nMelBins), Int32(nFreqs), 1.0, melFilterbank, Int32(nFreqs), |
| 141 | + powerSpec, 1, 0.0, &melFrame, 1) |
| 142 | + |
| 143 | + for i in 0..<nMelBins { |
| 144 | + mel[i * nFrames + t] = log10(max(melFrame[i], 1e-10)) |
| 145 | + } |
| 146 | + } |
| 147 | + |
| 148 | + // Normalise: clamp to max−8, then (x+4)/4 |
| 149 | + let maxVal = mel.max() ?? 0 |
| 150 | + for i in 0..<mel.count { mel[i] = (max(mel[i], maxVal - 8) + 4) / 4 } |
| 151 | + return mel |
| 152 | + } |
| 153 | + |
| 154 | + // MARK: - Filterbank builder |
| 155 | + |
| 156 | + private static func makeMelFilterbank() -> [Float] { |
| 157 | + let fMax: Float = Float(sampleRate) / 2 // 8000 Hz |
| 158 | + |
| 159 | + func hzToMel(_ f: Float) -> Float { 2595 * log10(1 + f / 700) } |
| 160 | + func melToHz(_ m: Float) -> Float { 700 * (pow(10, m / 2595) - 1) } |
| 161 | + |
| 162 | + let melMin = hzToMel(0), melMax = hzToMel(fMax) |
| 163 | + let nPts = nMelBins + 2 |
| 164 | + let pts = (0..<nPts).map { i -> Float in |
| 165 | + melToHz(melMin + Float(i) / Float(nPts - 1) * (melMax - melMin)) |
| 166 | + } |
| 167 | + // FFT bin frequencies for n_fft = 400 |
| 168 | + let fftFreqs = (0..<nFreqs).map { Float($0) * Float(sampleRate) / Float(nFFT) } |
| 169 | + |
| 170 | + var fb = [Float](repeating: 0, count: nMelBins * nFreqs) |
| 171 | + for m in 0..<nMelBins { |
| 172 | + let fL = pts[m], fC = pts[m + 1], fR = pts[m + 2] |
| 173 | + let norm: Float = 2 / (fR - fL) |
| 174 | + for k in 0..<nFreqs { |
| 175 | + let f = fftFreqs[k] |
| 176 | + if f >= fL && f <= fC { fb[m * nFreqs + k] = norm * (f - fL) / (fC - fL) } |
| 177 | + else if f > fC && f <= fR { fb[m * nFreqs + k] = norm * (fR - f) / (fR - fC) } |
| 178 | + } |
| 179 | + } |
| 180 | + return fb |
| 181 | + } |
| 182 | +} |
0 commit comments