synthesis bugfix

Project.toml

@@ -1,7 +1,7 @@
 name = "STFT"
 uuid = "58bb99bf-048b-48b7-93e7-1cbf3ee61509"
 authors = ["Szymon M. Woźniak <s@zymon.org>"]
-version = "1.2.0"
+version = "1.3.0"
 
 [deps]
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"

src/STFT.jl

@@ -5,9 +5,21 @@ using FFTW
 
 export stft, istft
 
-_fft(x::AbstractMatrix{<:Real}, d) = rfft(x, d)
+_fft(x::AbstractArray{<:Real}, d) = rfft(x, d)
-_fft(x::AbstractMatrix{<:Complex}, d) = fft(x, d)
+_fft(x::AbstractArray{<:Complex}, d) = fft(x, d)
 
+function malloc_stft(
+    x::M,
+    w::V,
+    L::I = zero(I),
+    N::I = length(w);
+) where {T<:Number, I<:Integer, V<:AbstractVector{T}, M<:AbstractMatrix{T}}
+    X, K = size(x)      # Length of the signal in samples
+    W = length(w)       # Length of the window in samples
+    S = (X-L) ÷ (W - L) # Number of segments
+    N = N < W ? W : N   # DFT size
+    zeros(T, (N, S, K)) # Allocate container for signal segments
+end
 
 
 doc_analysis = """
@@ -94,34 +106,32 @@ function analysis() end
 "$doc_analysis"
 stft(x, w, L=0, N=length(w)) = analysis(x, w, L, N)
 
-function analysis(
-    x::V,
-    w::V,
-    L::I = zero(I),
-    N::I = length(w);
-)::Matrix{T |> complex} where {T<:Number, I<:Integer, V<:AbstractVector{T}}
-    analysis((@view x[:,:]), w, L, N)
-end
-
 function analysis(
     x::M,
     w::V,
     L::I = zero(I),
     N::I = length(w);
 ) where {T<:Number, I<:Integer, V<:AbstractVector{T}, M<:AbstractMatrix{T}}
-    X, K = size(x)      # Length of the signal in samples
+    sc = malloc_stft(x, w, L, N)
-    W = length(w)       # Length of the window in samples
+    N, S, K = sc |> size
-    H = W - L           # Hop
+    W = w |> length
-    S = (X-L) ÷ H       # Number of segments
-    N = N < W ? W : N   # DFT size
-    sc = zeros(T, N, S, K) # Allocate container for signal segments
 
-    @turbo for s ∈ 1:S, k ∈ 1:K, n ∈ 1:W
+    @turbo for k ∈ 1:K, s ∈ 1:S, n ∈ 1:W
-        sc[n, s, k] = w[n] * x[(s-1)*H+n, k]
+        sc[n, s, k] = w[n] * x[(s-1)*(W-L)+n, k]
     end
     _fft(sc, 1) # Convert segments to frequency-domain
 end
 
+function analysis(
+    x::V,
+    w::V,
+    L::I = zero(I),
+    N::I = length(w);
+)::Matrix{T |> complex} where {T<:Number, I<:Integer, V<:AbstractVector{T}}
+    xx = @view x[:,:]
+    @view analysis(xx, w, L, N)[:, :, 1]
+end
+
 function analysis(
     xs::AbstractArray{<:AbstractVector{T}},
     w::AbstractVector{T},
@@ -239,14 +249,32 @@ function synthesis(
 
     xs = irfft(X, N, 1) # Convert segments to time-domain
 
-    @turbo for s ∈ 1:S
+    @turbo for s ∈ 0:(S-1)
-        ss = (s-1)*H # Segment start
+        ss = s*H # Segment shift
         for n = 1:W
-            xn[ss+n] += xs[n, s] * w[n]
+            xn[ss+n] += xs[n, s+1] * w[n]
             xd[ss+n] += w²[n]
         end
     end
     xn ./ xd # Normalize
 end
 
+"""
+Real-value signal STFT with constant input size.
+"""
+function rSTFTm(A, w, L, N=length(w))
+    mem = STFT.malloc_stft(A, w, L, N)
+    N, S, K = mem |> size
+    W = w |> length
+    P = plan_rfft(mem, 1)
+    function f(x)
+        @turbo for k ∈ 1:K, s ∈ 1:S, n ∈ 1:W
+            mem[n, s, k] = w[n] * x[(s-1)*(W-L)+n, k]
+        end
+        P * mem
+    end
+    return f
+end
+
+
 end # module