0001    // Convolution.swift
0002    //
0003    // Copyright (c) 2014–2015 Mattt Thompson (http://mattt.me)
0004    // Copyright (c) 2015-2016 Remy Prechelt
0005    //
0006    // Permission is hereby granted, free of charge, to any person obtaining a copy
0007    // of this software and associated documentation files (the "Software"), to deal
0008    // in the Software without restriction, including without limitation the rights
0009    // to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
0010    // copies of the Software, and to permit persons to whom the Software is
0011    // furnished to do so, subject to the following conditions:
0012    //
0013    // The above copyright notice and this permission notice shall be included in
0014    // all copies or substantial portions of the Software.
0015    //
0016    // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
0017    // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
0018    // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
0019    // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
0020    // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
0021    // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
0022    // THE SOFTWARE.
0023    
0024    import Accelerate
0025    
0026    // MARK: Convolution
0027    
0028    // Convolution of a signal [x], with a kernel [k]. The signal must be at least as long as the kernel.
0029    public func conv(x: [Float], _ k: [Float]) -> [Float] {
0030        precondition(x.count >= k.count, "Input vector [x] must have at least as many elements as the kernel,  [k]")
0031        
0032        let resultSize = x.count + k.count - 1
0033        var result = [Float](count: resultSize, repeatedValue: 0)
0034        let kEnd = UnsafePointer<Float>(k).advancedBy(k.count - 1)
0035        let xPad = Repeat(count: k.count-1, repeatedValue: Float(0.0))
0036        let xPadded = xPad + x + xPad
0037        vDSP_conv(xPadded, 1, kEnd, -1, &result, 1, vDSP_Length(resultSize), vDSP_Length(k.count))
0038    
0039        return result
0040    }
0041    
0042    // Convolution of a signal [x], with a kernel [k]. The signal must be at least as long as the kernel.
0043    public func conv(x: [Double], _ k: [Double]) -> [Double] {
0044        precondition(x.count >= k.count, "Input vector [x] must have at least as many elements as the kernel,  [k]")
0045        
0046        let resultSize = x.count + k.count - 1
0047        var result = [Double](count: resultSize, repeatedValue: 0)
0048        let kEnd = UnsafePointer<Double>(k).advancedBy(k.count - 1)
0049        let xPad = Repeat(count: k.count-1, repeatedValue: Double(0.0))
0050        let xPadded = xPad + x + xPad
0051        vDSP_convD(xPadded, 1, kEnd, -1, &result, 1, vDSP_Length(resultSize), vDSP_Length(k.count))
0052        
0053        return result
0054    }
0055    
0056    // MARK: Cross-Correlation
0057    
0058    // Cross-correlation of a signal [x], with another signal [y]. The signal [y]
0059    // is padded so that it is the same length as [x].
0060    public func xcorr(x: [Float], _ y: [Float]) -> [Float] {
0061        precondition(x.count >= y.count, "Input vector [x] must have at least as many elements as [y]")
0062        var yPadded = y
0063        if x.count > y.count {
0064            let padding = Repeat(count: x.count - y.count, repeatedValue: Float(0.0))
0065            yPadded = y + padding
0066        }
0067        
0068        let resultSize = x.count + yPadded.count - 1
0069        var result = [Float](count: resultSize, repeatedValue: 0)
0070        let xPad = Repeat(count: yPadded.count-1, repeatedValue: Float(0.0))
0071        let xPadded = xPad + x + xPad
0072        vDSP_conv(xPadded, 1, yPadded, 1, &result, 1, vDSP_Length(resultSize), vDSP_Length(yPadded.count))
0073        
0074        return result
0075    }
0076    
0077    // Cross-correlation of a signal [x], with another signal [y]. The signal [y]
0078    // is padded so that it is the same length as [x].
0079    public func xcorr(x: [Double], _ y: [Double]) -> [Double] {
0080        precondition(x.count >= y.count, "Input vector [x] must have at least as many elements as [y]")
0081        var yPadded = y
0082        if x.count > y.count {
0083            let padding = Repeat(count: x.count - y.count, repeatedValue: Double(0.0))
0084            yPadded = y + padding
0085        }
0086        
0087        let resultSize = x.count + yPadded.count - 1
0088        var result = [Double](count: resultSize, repeatedValue: 0)
0089        let xPad = Repeat(count: yPadded.count-1, repeatedValue: Double(0.0))
0090        let xPadded = xPad + x + xPad
0091        vDSP_convD(xPadded, 1, yPadded, 1, &result, 1, vDSP_Length(resultSize), vDSP_Length(yPadded.count))
0092        
0093        return result
0094    }
0095    
0096    // MARK: Auto-correlation
0097    
0098    // Auto-correlation of a signal [x]
0099    public func xcorr(x: [Float]) -> [Float] {
0100        let resultSize = 2*x.count - 1
0101        var result = [Float](count: resultSize, repeatedValue: 0)
0102        let xPad = Repeat(count: x.count-1, repeatedValue: Float(0.0))
0103        let xPadded = xPad + x + xPad
0104        vDSP_conv(xPadded, 1, x, 1, &result, 1, vDSP_Length(resultSize), vDSP_Length(x.count))
0105        
0106        return result
0107    }
0108    
0109    // Auto-correlation of a signal [x]
0110    public func xcorr(x: [Double]) -> [Double] {
0111        let resultSize = 2*x.count - 1
0112        var result = [Double](count: resultSize, repeatedValue: 0)
0113        let xPad = Repeat(count: x.count-1, repeatedValue: Double(0.0))
0114        let xPadded = xPad + x + xPad
0115        vDSP_convD(xPadded, 1, x, 1, &result, 1, vDSP_Length(resultSize), vDSP_Length(x.count))
0116        
0117        return result
0118    }
0119