Skip to content

Fourier, DCT and Hartley Transforms

Fast transforms on real data, including zoom FFT and FFT-based convolution.

Fourier

v_rfft

V_RFFT - Calculate the DFT of real data, returning only the first half.

v_rfft 

v_rfft(x, n=None, d=None) -> ndarray

Calculate the DFT of real data Y=(X,N,D).

Data is truncated/padded to length N if specified. N even: (N+2)/2 points are returned with the first and last being real N odd: (N+1)/2 points are returned with the first being real In all cases floor(1+N/2) points are returned. D is the dimension (0-based axis) along which to do the DFT.

Parameters:

Name Type Description Default
x array_like

Real input data.

 required
n int

Transform length. Default is the size of x along axis d.

 None
d int

Axis along which to compute the FFT. Default is first axis with size > 1.

 None

Returns:

Name Type Description
y ndarray

The first floor(1+N/2) points of the DFT.
Source code in pyvoicebox/v_rfft.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
def v_rfft(x, n=None, d=None) -> np.ndarray:
    """Calculate the DFT of real data Y=(X,N,D).

    Data is truncated/padded to length N if specified.
      N even: (N+2)/2 points are returned with the first and last being real
      N odd:  (N+1)/2 points are returned with the first being real
    In all cases floor(1+N/2) points are returned.
    D is the dimension (0-based axis) along which to do the DFT.

    Parameters
    ----------
    x : array_like
        Real input data.
    n : int, optional
        Transform length. Default is the size of x along axis d.
    d : int, optional
        Axis along which to compute the FFT. Default is first axis with size > 1.

    Returns
    -------
    y : ndarray
        The first floor(1+N/2) points of the DFT.
    """
    x = np.asarray(x, dtype=float)
    s = x.shape

    # Scalar case
    if x.size == 1:
        return x.copy()

    if d is None:
        # Find first non-singleton dimension
        for i, si in enumerate(s):
            if si > 1:
                d = i
                break
        if d is None:
            d = 0
        if n is None:
            n = s[d]

    if n is None:
        n = s[d]

    y = np.fft.fft(x, n=n, axis=d)

    # Keep only first floor(1+n/2) points along axis d
    keep = 1 + n // 2
    slices = [slice(None)] * y.ndim
    slices[d] = slice(0, keep)
    y = y[tuple(slices)]

    return y

v_irfft

V_IRFFT - Inverse FFT of a conjugate symmetric spectrum.

v_irfft 

v_irfft(y, n=None, d=None) -> ndarray

Inverse FFT of a conjugate symmetric spectrum X=(Y,N,D).

Parameters:

Name Type Description Default
y array_like

The first half of a complex spectrum (as produced by v_rfft).

 required
n int

Number of output points to generate. Default: 2*M-2 where M is size of y along axis d. Note: default is always even; specify n explicitly if odd output is desired.

 None
d int

Axis along which to perform the transform. Default: first non-singleton dimension.

 None

Returns:

Name Type Description
x ndarray

Real inverse DFT of y.
Source code in pyvoicebox/v_irfft.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
def v_irfft(y, n=None, d=None) -> np.ndarray:
    """Inverse FFT of a conjugate symmetric spectrum X=(Y,N,D).

    Parameters
    ----------
    y : array_like
        The first half of a complex spectrum (as produced by v_rfft).
    n : int, optional
        Number of output points to generate. Default: 2*M-2 where M is size
        of y along axis d. Note: default is always even; specify n explicitly
        if odd output is desired.
    d : int, optional
        Axis along which to perform the transform. Default: first non-singleton
        dimension.

    Returns
    -------
    x : ndarray
        Real inverse DFT of y.
    """
    y = np.asarray(y, dtype=complex)
    s = list(y.shape)
    ns = len(s)

    # Scalar case
    if y.size == 1:
        return np.real(y).copy()

    if d is None:
        for i, si in enumerate(s):
            if si > 1:
                d = i
                break
        if d is None:
            d = 0

    m = s[d]
    k = y.size // m  # number of FFTs to do

    # Reshape: move dimension d to the front, then reshape to (m, k)
    if d == 0:
        v = y.reshape(m, k)
    else:
        perm = list(range(d, ns)) + list(range(0, d))
        v = np.transpose(y, perm).reshape(m, k)

    if n is None:
        n = 2 * m - 2  # default output length
    else:
        mm = 1 + n // 2  # expected input length
        if mm > m:
            v = np.vstack([v, np.zeros((mm - m, k), dtype=complex)])
        elif mm < m:
            v = v[:mm, :]
        m = mm

    v[0, :] = np.real(v[0, :])  # force DC element real

    if n % 2:  # odd output length
        full = np.vstack([v, np.conj(v[m - 1:0:-1, :])])
        x = np.real(np.fft.ifft(full, axis=0))
    else:  # even output length
        v[m - 1, :] = np.real(v[m - 1, :])  # force Nyquist element real
        w = np.ones((1, k))
        t = -0.5j * np.exp((2j * np.pi / n) * np.arange(m))
        t = t[:, np.newaxis]
        z = (t + 0.5) * (np.conj(v[::-1, :]) - v) + v
        z = z[:m - 1, :]  # remove last row (z[m-1,:])
        zz = np.fft.ifft(z, axis=0)
        x = np.zeros((n, k))
        x[0::2, :] = np.real(zz)
        x[1::2, :] = np.imag(zz)

    s[d] = n  # change output dimension
    if d == 0:
        x = x.reshape(s)
    else:
        # Reorder dimensions: s was permuted as [d:ns, 0:d]
        perm_shape = [s[i] for i in (list(range(d, ns)) + list(range(0, d)))]
        x = x.reshape(perm_shape)
        # Inverse permutation
        inv_perm = list(range(ns + 1 - d, ns)) + list(range(0, ns + 1 - d))
        x = np.transpose(x, inv_perm)

    return x

v_rsfft

V_RSFFT - FFT of a real symmetric spectrum.

v_rsfft 

v_rsfft(y, n=None) -> ndarray

FFT of a real symmetric spectrum X=(Y,N).

Y is the "first half" of a symmetric real input signal and X is the "first half" of the symmetric real Fourier transform.

If N is even, the "first half" contains 1+N/2 elements. If N is odd, the "first half" contains (N+1)/2 elements. If N is omitted it will be taken to be 2*(length(Y)-1) and is always even.

If Y is a matrix (2-D), the transform is performed along each column.

The inverse function is y = v_rsfft(x, n) / n.

Parameters:

Name Type Description Default
y array_like

Real input data. Must be real-valued.

 required
n int

Full signal length. Default: 2*(M-1) where M is the number of rows.

 None

Returns:

Name Type Description
x ndarray

The first half of the symmetric real Fourier transform.
Source code in pyvoicebox/v_rsfft.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
def v_rsfft(y, n=None) -> np.ndarray:
    """FFT of a real symmetric spectrum X=(Y,N).

    Y is the "first half" of a symmetric real input signal and X is the
    "first half" of the symmetric real Fourier transform.

    If N is even, the "first half" contains 1+N/2 elements.
    If N is odd, the "first half" contains (N+1)/2 elements.
    If N is omitted it will be taken to be 2*(length(Y)-1) and is always even.

    If Y is a matrix (2-D), the transform is performed along each column.

    The inverse function is y = v_rsfft(x, n) / n.

    Parameters
    ----------
    y : array_like
        Real input data. Must be real-valued.
    n : int, optional
        Full signal length. Default: 2*(M-1) where M is the number of rows.

    Returns
    -------
    x : ndarray
        The first half of the symmetric real Fourier transform.
    """
    y = np.asarray(y, dtype=float)
    if not np.all(np.isreal(y)):
        raise ValueError('RSFFT: Input must be real')

    fl = False
    if y.ndim == 1:
        fl = True
        y = y[:, np.newaxis]
    elif y.shape[0] == 1:
        fl = True
        y = y.T

    m, k = y.shape

    if n is None:
        n = 2 * m - 2
    else:
        mm = 1 + n // 2
        if mm > m:
            y = np.vstack([y, np.zeros((mm - m, k))])
        elif mm < m:
            y = y[:mm, :]
        m = mm

    # Build full symmetric signal and take FFT
    # MATLAB: fft([y; y(n-m+1:-1:2, :)])
    # n-m+1 in 1-based -> n-m in 0-based, down to index 2 in 1-based -> index 1 in 0-based
    full = np.vstack([y, y[n - m:0:-1, :]])
    x = np.real(np.fft.fft(full, axis=0))
    x = x[:m, :]  # keep first m rows

    if fl:
        x = x.ravel()

    return x

v_zoomfft

V_ZOOMFFT - DTFT evaluated over a linear frequency range.

v_zoomfft 

v_zoomfft(
    x, n=None, m=None, s=0, d=None
) -> tuple[ndarray, ndarray]

DTFT evaluated over a linear frequency range Y=(X,N,M,S,D).

Parameters:

Name Type Description Default
x array_like

Input vector or array.

 required
n float

Reciprocal of normalized frequency increment (can be non-integer). The frequency increment is fs/n. Default: size(x, d).

 None
m int

Number of output points is floor(m). Default: floor(n).

 None
s float

Starting frequency index (can be non-integer). The starting frequency is s*fs/n. Default: 0.

 0
d int

Axis along which to do FFT (0-based). Default: first non-singleton.

 None

Returns:

Name Type Description
y ndarray

Output DTFT coefficients.
f ndarray

Normalized frequencies (1 corresponds to fs).
Source code in pyvoicebox/v_zoomfft.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
def v_zoomfft(x, n=None, m=None, s=0, d=None) -> tuple[np.ndarray, np.ndarray]:
    """DTFT evaluated over a linear frequency range Y=(X,N,M,S,D).

    Parameters
    ----------
    x : array_like
        Input vector or array.
    n : float, optional
        Reciprocal of normalized frequency increment (can be non-integer).
        The frequency increment is fs/n. Default: size(x, d).
    m : int, optional
        Number of output points is floor(m). Default: floor(n).
    s : float, optional
        Starting frequency index (can be non-integer).
        The starting frequency is s*fs/n. Default: 0.
    d : int, optional
        Axis along which to do FFT (0-based). Default: first non-singleton.

    Returns
    -------
    y : ndarray
        Output DTFT coefficients.
    f : ndarray
        Normalized frequencies (1 corresponds to fs).
    """
    x = np.asarray(x, dtype=complex)
    e = list(x.shape)
    p = x.size

    if d is None:
        dims = [i for i, ei in enumerate(e) if ei > 1]
        if dims:
            d = dims[0]
        else:
            d = 0

    k_len = e[d]
    q = p // k_len

    if d == 0:
        z = x.reshape(k_len, q)
    else:
        z = np.moveaxis(x, d, 0).reshape(k_len, q)

    if n is None:
        n = k_len
    if m is None:
        m = int(np.floor(n))
    else:
        m = int(np.floor(m))
    if s is None:
        s = 0

    k = k_len

    l = int(2 ** np.ceil(np.log2(m + k - 1)))  # next power of 2

    if n == int(n) and s == int(s) and int(n) < 2 * l and int(n) >= k:
        # Quickest to do a normal FFT
        ni = int(n)
        a = np.fft.fft(z, n=ni, axis=0)
        si = int(s)
        indices = np.mod(np.arange(si, si + m), ni)
        y = a[indices, :]
    else:
        # Chirp z-transform (Bluestein's algorithm)
        b = np.exp(1j * np.pi * np.mod((s + np.arange(1 - k, m)) ** 2, 2 * n) / n)
        c = np.conj(b[k - 1:k - 1 + m])
        h = np.fft.fft(b, n=l)
        g = np.exp(-1j * np.pi * np.mod(np.arange(k) ** 2, 2 * n) / n)
        a = np.fft.ifft(np.fft.fft(z * g[:, np.newaxis], n=l, axis=0) * h[:, np.newaxis], axis=0)
        y = a[k - 1:k - 1 + m, :] * c[:, np.newaxis]

    if d == 0:
        e_out = list(e)
        e_out[d] = m
        y = y.reshape(e_out)
    else:
        moved_shape = list(np.moveaxis(x, d, 0).shape)
        moved_shape[0] = m
        y = y.reshape(moved_shape)
        y = np.moveaxis(y, 0, d)

    f = (s + np.arange(m)) / n

    return y, f

DCT and Hartley

v_rdct

V_RDCT - Discrete cosine transform of real data.

v_rdct 

v_rdct(x, n=None, a=None, b=1.0) -> ndarray

Discrete cosine transform of real data Y=(X,N,A,B).

Parameters:

Name Type Description Default
x array_like

Real-valued input data. Transform is applied to columns.

 required
n int

Transform length. x will be zero-padded or truncated to length n. Default: number of rows of x.

 None
a float

Output scale factor. Default: sqrt(2*n).

 None
b float

Output scale factor for DC term. Default: 1.

 1.0

Returns:

Name Type Description
y ndarray

DCT output data.
Source code in pyvoicebox/v_rdct.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
def v_rdct(x, n=None, a=None, b=1.0) -> np.ndarray:
    """Discrete cosine transform of real data Y=(X,N,A,B).

    Parameters
    ----------
    x : array_like
        Real-valued input data. Transform is applied to columns.
    n : int, optional
        Transform length. x will be zero-padded or truncated to length n.
        Default: number of rows of x.
    a : float, optional
        Output scale factor. Default: sqrt(2*n).
    b : float, optional
        Output scale factor for DC term. Default: 1.

    Returns
    -------
    y : ndarray
        DCT output data.
    """
    x = np.asarray(x, dtype=float)

    fl = False
    if x.ndim == 1:
        fl = True
        x = x[:, np.newaxis]
    elif x.shape[0] == 1:
        fl = True
        x = x.T

    m, k = x.shape
    if n is None:
        n = m
    if a is None:
        a = np.sqrt(2.0 * n)

    # Zero-pad or truncate
    if n > m:
        x = np.vstack([x, np.zeros((n - m, k))])
    elif n < m:
        x = x[:n, :]

    # Reorder: [x(1:2:n,:); x(2*fix(n/2):-2:2,:)]
    # MATLAB 1-based: odd indices 1,3,5,... then even indices from 2*fix(n/2) down by 2 to 2
    odd_idx = np.arange(0, n, 2)  # 0-based: 0,2,4,...
    even_top = 2 * (n // 2) - 1   # 0-based version of 2*fix(n/2)
    even_idx = np.arange(even_top, 0, -2)  # 0-based: ...,3,1
    idx = np.concatenate([odd_idx, even_idx])
    x = x[idx, :]

    # Compute z multipliers
    # MATLAB: z = [sqrt(2) 2*exp((-0.5i*pi/n)*(1:n-1))].';
    z = np.zeros(n, dtype=complex)
    z[0] = np.sqrt(2.0)
    z[1:] = 2.0 * np.exp((-0.5j * np.pi / n) * np.arange(1, n))

    # y = real(fft(x) .* z(:, ones(1,k))) / a
    y = np.real(np.fft.fft(x, axis=0) * z[:, np.newaxis]) / a
    y[0, :] = y[0, :] * b

    if fl:
        y = y.ravel()

    return y

v_irdct

V_IRDCT - Inverse discrete cosine transform of real data.

v_irdct 

v_irdct(y, n=None, a=None, b=1.0) -> ndarray

Inverse discrete cosine transform of real data X=(Y,N,A,B).

Parameters:

Name Type Description Default
y array_like

DCT coefficients (as produced by v_rdct). Transform applied to columns.

 required
n int

Output length. Default: number of rows of y.

 None
a float

Scale factor. Default: sqrt(2*m) where m is the number of rows of y.

 None
b float

DC scale factor. Default: 1.

 1.0

Returns:

Name Type Description
x ndarray

Inverse DCT output data.
Source code in pyvoicebox/v_irdct.py 
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
def v_irdct(y, n=None, a=None, b=1.0) -> np.ndarray:
    """Inverse discrete cosine transform of real data X=(Y,N,A,B).

    Parameters
    ----------
    y : array_like
        DCT coefficients (as produced by v_rdct). Transform applied to columns.
    n : int, optional
        Output length. Default: number of rows of y.
    a : float, optional
        Scale factor. Default: sqrt(2*m) where m is the number of rows of y.
    b : float, optional
        DC scale factor. Default: 1.

    Returns
    -------
    x : ndarray
        Inverse DCT output data.
    """
    y = np.asarray(y, dtype=float)

    fl = False
    if y.ndim == 1:
        fl = True
        y = y[:, np.newaxis]
    elif y.shape[0] == 1:
        fl = True
        y = y.T

    m_orig, k = y.shape

    if a is None:
        a = np.sqrt(2.0 * m_orig)
    if n is None:
        n = m_orig

    # Zero-pad or truncate
    if n > m_orig:
        y = np.vstack([y, np.zeros((n - m_orig, k))])
    elif n < m_orig:
        y = y[:n, :]

    x = np.zeros((n, k))
    m = (n + 1) // 2  # fix((n+1)/2)
    p = n - m

    # z = 0.5*exp((0.5i*pi/n)*(1:p)).'
    z = 0.5 * np.exp((0.5j * np.pi / n) * np.arange(1, p + 1))
    z = z[:, np.newaxis]

    # u = (y(2:p+1,:) - 1i*y(n:-1:m+1,:)) .* z(:,w) * a
    # MATLAB 1-based: y(2:p+1,:) -> 0-based: y[1:p+1,:]
    # MATLAB 1-based: y(n:-1:m+1,:) -> 0-based: y[n-1:m:-1,:]  (but n is already length, so y has indices 0..n-1)
    u = (y[1:p + 1, :] - 1j * y[n - 1:m - 1:-1, :]) * z * a

    # y = [y(1,:)*sqrt(0.5)*a/b; u(1:m-1,:)]
    yy = np.vstack([y[0:1, :] * np.sqrt(0.5) * a / b, u[:m - 1, :]])

    if m == p:
        # Even n case
        # z = -0.5i*exp((2i*pi/n)*(0:m-1)).'
        z2 = -0.5j * np.exp((2j * np.pi / n) * np.arange(m))
        z2 = z2[:, np.newaxis]
        # y = (z(:,w)+0.5).*(conj(flipud(u))-y)+y
        yy = (z2 + 0.5) * (np.conj(u[::-1, :]) - yy) + yy
        # z = ifft(y,[],1)
        zz = np.fft.ifft(yy, axis=0)
        uu = np.real(zz)
        yi = np.imag(zz)
        q = m // 2
        h = (m % 2) / 2.0  # rem(m,2)/2

        # MATLAB uses 1-based indexing for x(1:4:n,:), x(2:4:n,:), etc.
        # In 0-based: x[0::4], x[1::4], x[2::4], x[3::4]
        # MATLAB: x(1:4:n,:)=u(1:q+h,:)   -> 0-based: x[0::4]=uu[0:int(q+h)]
        # MATLAB: x(2:4:n,:)=y(m:-1:q+1-h,:) -> 0-based: x[1::4]=yi[m-1:int(q-h):-1]
        # MATLAB: x(3:4:n,:)=y(1:q-h,:) -> 0-based: x[2::4]=yi[0:int(q-h)]
        # MATLAB: x(4:4:n,:)=u(m:-1:q+1+h,:) -> 0-based: x[3::4]=uu[m-1:int(q+h):-1]
        # Note: h is 0 when m is even, 0.5 when m is odd
        # But these are used as integer indices; in MATLAB q+h is floor/ceil depending
        ih = int(h)  # 0 if m even, but h=0.5 rounds... let's be careful

        # When m is even: h=0, q=m/2
        #   x(1:4:n)=u(1:q) -> x[0::4]=uu[0:q]
        #   x(2:4:n)=y(m:-1:q+1) -> x[1::4]=yi[m-1:q-1:-1]
        #   x(3:4:n)=y(1:q) -> x[2::4]=yi[0:q]
        #   x(4:4:n)=u(m:-1:q+1) -> x[3::4]=uu[m-1:q-1:-1]
        # When m is odd: h=0.5, q=(m-1)/2
        #   x(1:4:n)=u(1:q+0.5)=u(1:q+1) (since MATLAB rounds up for indexing)
        #   Actually in MATLAB, q+h where h=0.5 gives q+0.5, and 1:q+0.5 means 1:floor(q+0.5)
        #   For m odd, q=(m-1)/2, so q+0.5 = m/2 which is not integer... MATLAB truncates
        #   Actually MATLAB: 1:q+h where q=int, h=0.5 gives 1:q (since q+0.5 is not >= q+1)
        #   Wait, let me re-examine. In MATLAB, 1:2.5 = [1, 2], so it goes up to floor(2.5)=2.
        #   So: when m is odd, h=0.5:
        #     q = (m-1)/2  (integer since m is odd)
        #     q+h = q + 0.5 = m/2 (non-integer)
        #     1:q+h = 1:q  (MATLAB truncates to integer steps)
        #     q+1-h = q+0.5 (non-integer)
        #     m:-1:q+1-h = m:-1:q+1 (MATLAB ceil for reverse: goes down to ceil(q+0.5)=q+1)
        #     q-h = q-0.5 (non-integer)
        #     1:q-h = 1:q-1 (MATLAB floor)
        #     q+1+h = q+1.5
        #     m:-1:q+1+h = m:-1:q+2 (MATLAB ceil for reverse: ceil(q+1.5)=q+2)

        if m % 2 == 0:
            # m even
            x[0::4, :] = uu[:q, :]
            x[1::4, :] = yi[m - 1:q - 1:-1, :]
            x[2::4, :] = yi[:q, :]
            x[3::4, :] = uu[m - 1:q - 1:-1, :]
        else:
            # m odd: q = (m-1)//2
            x[0::4, :] = uu[:q, :]            # 1:q in MATLAB (q elements)
            x[1::4, :] = yi[m - 1:q:-1, :]    # m:-1:q+1 in MATLAB
            x[2::4, :] = yi[:q, :]            # 1:q-1 -> but wait, need to check count
            # Actually let me count more carefully for odd m
            # n = 2*m when m==p, so n is even
            # x has n elements, x[0::4] has n/4 elements = m/2 elements
            # For m odd, that's (m-1)/2 = q... but len(x[0::4]) could be ceil(n/4)
            # n = 2m, so x[0::4] has ceil(2m/4) = ceil(m/2) elements
            # For m odd, ceil(m/2) = (m+1)/2 = q+1
            x[0::4, :] = uu[:q + 1, :]          # q+1 elements
            x[1::4, :] = yi[m - 1:q:-1, :]      # m-1 down to q+1 = m-1-q elements
            x[2::4, :] = yi[:q, :]              # q elements
            x[3::4, :] = uu[m - 1:q:-1, :]      # m-1 down to q+1 = m-1-q elements
    else:
        # Odd n case
        # z = real(ifft([y; conj(flipud(u))]))
        full = np.vstack([yy, np.conj(u[::-1, :])])
        zz = np.real(np.fft.ifft(full, axis=0))
        x[0::2, :] = zz[:m, :]
        x[1::2, :] = zz[n - 1:m - 1:-1, :]

    if fl:
        x = x.ravel()

    return x

v_rhartley

V_RHARTLEY - Calculate the Hartley transform of real data.

v_rhartley 

v_rhartley(x, n=None) -> ndarray

Calculate the Hartley transform of real data Y=(X,N).

Data is truncated/padded to length N if specified. The inverse transformation is x = v_rhartley(y, n) / n.

Parameters:

Name Type Description Default
x array_like

Real input data.

 required
n int

Transform length. Default: length of x.

 None

Returns:

Name Type Description
y ndarray

Hartley transform of x.
Source code in pyvoicebox/v_rhartley.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
def v_rhartley(x, n=None) -> np.ndarray:
    """Calculate the Hartley transform of real data Y=(X,N).

    Data is truncated/padded to length N if specified.
    The inverse transformation is x = v_rhartley(y, n) / n.

    Parameters
    ----------
    x : array_like
        Real input data.
    n : int, optional
        Transform length. Default: length of x.

    Returns
    -------
    y : ndarray
        Hartley transform of x.
    """
    x = np.asarray(x, dtype=float)
    x = np.real(x)

    # MATLAB's fft operates on the first non-singleton dimension by default
    # (columns for 2-D arrays). Match this behavior.
    if x.ndim >= 2:
        axis = 0
    else:
        axis = -1

    if n is None:
        y = np.fft.fft(x, axis=axis)
    else:
        y = np.fft.fft(x, n=n, axis=axis)

    y = np.real(y) - np.imag(y)

    return y

Convolution

v_convfft

V_CONVFFT - 1-D convolution or correlation using FFT.

v_convfft 

v_convfft(
    x, h, d=None, m="", h0=1, x1=1, x2=None
) -> ndarray

1-D convolution or correlation using FFT.

Parameters:

Name Type Description Default
x array_like or int

Input array, or size(x, d) when using precompute mode ('z' in m). If h is a _ConvFFTPrecomputed object, x is the input array.

 required
h array_like or _ConvFFTPrecomputed

Impulse response, or precomputed structure from previous call with 'z'.

 required
d int

Axis of x to do convolution along (0-based). Default: first non-singleton.

 None
m str

Mode options: 'x' for real correlation, 'X' for complex correlation, 'z' for precomputation.

 ''
h0 int

Origin sample number in h (1-based). Default: 1.

 1
x1 int

Start of range in x to align with origin of h (1-based). Default: 1.

 1
x2 int

End of range in x to align with origin of h (1-based). Default: size(x, d).

 None

Returns:

Name Type Description
z ndarray or _ConvFFTPrecomputed

Convolution output, or precomputed structure if 'z' in m.
Source code in pyvoicebox/v_convfft.py 
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
def v_convfft(x, h, d=None, m='', h0=1, x1=1, x2=None) -> np.ndarray:
    """1-D convolution or correlation using FFT.

    Parameters
    ----------
    x : array_like or int
        Input array, or size(x, d) when using precompute mode ('z' in m).
        If h is a _ConvFFTPrecomputed object, x is the input array.
    h : array_like or _ConvFFTPrecomputed
        Impulse response, or precomputed structure from previous call with 'z'.
    d : int, optional
        Axis of x to do convolution along (0-based). Default: first non-singleton.
    m : str, optional
        Mode options: 'x' for real correlation, 'X' for complex correlation,
        'z' for precomputation.
    h0 : int, optional
        Origin sample number in h (1-based). Default: 1.
    x1 : int, optional
        Start of range in x to align with origin of h (1-based). Default: 1.
    x2 : int, optional
        End of range in x to align with origin of h (1-based). Default: size(x, d).

    Returns
    -------
    z : ndarray or _ConvFFTPrecomputed
        Convolution output, or precomputed structure if 'z' in m.
    """
    if isinstance(h, _ConvFFTPrecomputed):
        # h is a precomputed structure
        precomp = h
        d = precomp.d
        nx = precomp.nx
        ns = precomp.ns
        nv = precomp.nv
        vmin = precomp.vmin
        vmax = precomp.vmax
        nf = precomp.nf
        fh = precomp.fh
        fmin = precomp.fmin
        fmax = precomp.fmax
        nz = precomp.nz
        zmin = precomp.zmin
        zmax = precomp.zmax

        x = np.asarray(x, dtype=complex)
        s = list(x.shape)
        k = x.size // nx
        if x.shape[d] != nx or len(s) != ns:
            raise ValueError('input x has incompatible dimensions')

        # Reshape x
        if d == 0:
            v = x.reshape(nx, k)
        else:
            perm = list(range(d, ns)) + list(range(0, d))
            v = np.transpose(x, perm).reshape(nx, k)

        if nv < nx:
            v = v[vmin:vmax + 1, :]  # 0-based inclusive

        # Do the convolution
        zz = np.fft.ifft(np.fft.fft(v, n=nf, axis=0) * fh[:, np.newaxis], axis=0)
        z_out = np.zeros((nz, k), dtype=complex)
        z_out[zmin:zmax + 1, :] = zz[fmin:fmax + 1, :]

        if np.all(np.isreal(x)):
            z_out = np.real(z_out)

        s[d] = nz
        if d == 0:
            z_out = z_out.reshape(s)
        else:
            perm_shape = [s[i] for i in (list(range(d, ns)) + list(range(0, d)))]
            z_out = z_out.reshape(perm_shape)
            inv_perm = list(range(ns - d, ns)) + list(range(0, ns - d))
            z_out = np.transpose(z_out, inv_perm)

        return z_out

    # Normal input calling sequence
    x = np.asarray(x)
    h = np.asarray(h, dtype=complex).ravel()

    s = list(x.shape)
    ps = x.size
    ns = len(s)

    if 'z' in m:
        # Output pre-computed structure
        if d is None:
            raise ValueError('d must be specified explicitly')
        if x.size != 1:
            raise ValueError('x must equal size(*, d)')
        nx = int(x.ravel()[0])
    else:
        if d is None:
            if ps < 2:
                d = 0
            else:
                for i, si in enumerate(s):
                    if si > 1:
                        d = i
                        break
                if d is None:
                    d = 0
        nx = s[d]

    k = ps // nx if ps > 1 else 1

    if x2 is None:
        x2 = nx

    nz = x2 - x1 + 1  # number of output lags
    nh = len(h)

    if 'X' in m:
        h = np.conj(h[::-1])
        h0 = nh + 1 - h0
    elif 'x' in m:
        h = h[::-1].copy()
        h0 = nh + 1 - h0

    hmin = h0 + x1 - nx
    hmax = h0 + x2 - 1
    xmin = x1 + h0 - nh
    xmax = x2 + h0 - 1

    if hmin > 1 or hmax < nh:
        hmin_c = max(hmin, 1)
        hmax_c = min(hmax, nh)
        h = h[hmin_c - 1:hmax_c]  # 1-based to 0-based
        nh = len(h)
        h0 = h0 - hmin_c + 1

    if xmin > 1 or xmax < nx:
        vmin = max(xmin, 1)
        vmax = min(xmax, nx)
        x1_new = x1 - vmin + h0
        x2_new = x2 - vmin + h0
    else:
        vmin = 1
        vmax = nx
        x1_new = x1 + h0 - 1
        x2_new = x2 + h0 - 1

    nv = vmax - vmin + 1
    nxz = min(max(max(nh - x1_new, 0), max(x2_new - nv, 0)), nh - 1)

    # Round up to next power of 2
    nf = int(2 ** np.ceil(np.log2(nv + nxz)))

    fmin_idx = max(x1_new, 1)
    fmax_idx = min(x2_new, min(nf, nx + nh - 1))
    zmin_idx = max(1, 2 - x1_new)
    zmax_idx = zmin_idx + fmax_idx - fmin_idx

    fh = np.fft.fft(h, n=nf)

    if 'z' in m:
        # Save as precomputed structure
        result = _ConvFFTPrecomputed()
        result.d = d
        result.nx = nx
        result.ns = ns
        result.nv = nv
        result.vmin = vmin - 1  # convert to 0-based
        result.vmax = vmax - 1  # convert to 0-based
        result.nf = nf
        result.fh = fh
        result.fmin = fmin_idx - 1  # convert to 0-based
        result.fmax = fmax_idx - 1  # convert to 0-based
        result.nz = nz
        result.zmin = zmin_idx - 1  # convert to 0-based
        result.zmax = zmax_idx - 1  # convert to 0-based
        return result

    if x.size > 0:
        # Reshape x
        if d == 0:
            v = x.reshape(nx, k)
        else:
            perm = list(range(d, ns)) + list(range(0, d))
            v = np.transpose(x, perm).reshape(nx, k)

        if nv < nx:
            v = v[vmin - 1:vmax, :]  # 1-based to 0-based

        v = np.asarray(v, dtype=complex)
        zz = np.fft.ifft(np.fft.fft(v, n=nf, axis=0) * fh[:, np.newaxis], axis=0)
        z_out = np.zeros((nz, k), dtype=complex)
        z_out[zmin_idx - 1:zmax_idx, :] = zz[fmin_idx - 1:fmax_idx, :]

        if np.all(np.isreal(x)):
            z_out = np.real(z_out)

        s[d] = nz
        if d == 0:
            z_out = z_out.reshape(s)
        else:
            perm_shape = [s[i] for i in (list(range(d, ns)) + list(range(0, d)))]
            z_out = z_out.reshape(perm_shape)
            inv_perm = list(range(ns - d, ns)) + list(range(0, ns - d))
            z_out = np.transpose(z_out, inv_perm)

        return z_out
    else:
        return np.array([])

v_frac2bin

V_FRAC2BIN - Convert a column vector to binary string representation.

v_frac2bin 

v_frac2bin(d, n=1, m=0) -> str

Convert a column vector to binary S=(D,N,M).

Parameters:

Name Type Description Default
d array_like

Scalar or 1-D array of values to convert.

 required
n int

Minimum number of integer bits to output. If negative, leading zeros will be output as spaces for positions to the left of |n|'th integer column. Default: 1.

 1
m int

Number of places after binary point. If negative, values are truncated rather than rounded. Default: 0.

 0

Returns:

Name Type Description
s list of str

List of binary string representations, one per input value.
Source code in pyvoicebox/v_frac2bin.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
def v_frac2bin(d, n=1, m=0) -> str:
    """Convert a column vector to binary S=(D,N,M).

    Parameters
    ----------
    d : array_like
        Scalar or 1-D array of values to convert.
    n : int, optional
        Minimum number of integer bits to output. If negative, leading zeros
        will be output as spaces for positions to the left of |n|'th integer
        column. Default: 1.
    m : int, optional
        Number of places after binary point. If negative, values are truncated
        rather than rounded. Default: 0.

    Returns
    -------
    s : list of str
        List of binary string representations, one per input value.
    """
    d = np.atleast_1d(np.asarray(d, dtype=float))
    l = abs(n)
    r = abs(m)

    # Find the maximum value's exponent
    max_val = np.max(d)
    if max_val > 0:
        _, e = np.frexp(max_val)  # e such that max_val = f * 2^e, 0.5 <= f < 1
    else:
        e = 0

    # Compute scaled values
    if m < 0:
        v = np.floor(d * (2.0 ** r))
    else:
        v = np.round(d * (2.0 ** r))

    v = v.astype(int)

    # Total number of bits needed: max(l, e) + r
    total_bits = max(l, e) + r

    # Generate binary strings
    result = []
    for val in v:
        if total_bits <= 0:
            bits = '0'
        else:
            bits = ''
            for bit_pos in range(total_bits - 1, -1, -1):
                bits += '1' if (val >> bit_pos) & 1 else '0'

        # Insert binary point if r > 0
        if r > 0:
            int_part = bits[:len(bits) - r]
            frac_part = bits[len(bits) - r:]
            if not int_part:
                int_part = '0'
            s = int_part + '.' + frac_part
        else:
            s = bits

        # Handle leading zeros -> spaces when n < 0
        if n < 0:
            bp = s.find('.')
            if bp < 0:
                bp = len(s)
            # Replace leading zeros with spaces for positions to the left of l'th column from point
            for i in range(bp - l):
                if s[i] == '0':
                    s = s[:i] + ' ' + s[i + 1:]
                else:
                    break

        result.append(s)

    return result