Быстрое косинусное преобразование через БПФ

15

Я хочу реализовать быстрое косинусное преобразование. Я прочитал в википедии , что есть быстрая версия DCT, которая аналогично вычисляется для FFT. Я попытался прочитать процитированную статью Makhoul * для реализаций FTPACK и FFTW, которые также используются в Scipy , но я не смог извлечь фактически алгоритм. Это то, что я до сих пор:

БПФ код:

def fft(x):
    if x.size ==1:
        return x
    N = x.size
    x0 = my_fft(x[0:N:2])
    x1 = my_fft(x[0+1:N:2])
    k = numpy.arange(N/2)
    e = numpy.exp(-2j*numpy.pi*k/N)
    l = x0 + x1 * e
    r = x0 - x1 * e  
    return numpy.hstack([l,r])

Код DCT:

def dct(x):
    k = 0
    N = x.size
    xk = numpy.zeros(N)
    for k in range(N):     
        for n in range(N):
            xn = x[n]
            xk[k] += xn*numpy.cos(numpy.pi/N*(n+1/2.0)*k)
    return xk 

Испытание FCT:

def my_fct(x):
    if x.size ==1:
        return x
    N = x.size
    x0 = my_fct(x[0:N:2]) # have to be set to zero?
    x1 = my_fct(x[0+1:N:2])
    k = numpy.arange(N/2)
    n = # ???
    c = numpy.cos(numpy.pi/N*(n+1/2.0)*k)
    l = x0 #???
    r = x0 #???
    return numpy.hstack([l,r])

* J. Махул, "Быстрое косинусное преобразование в одном и двух измерениях", IEEE Trans. Acoust. Речь Сиг. Proc. 28 (1), 27-34 (1980).

Framester
источник
2
Вы спрашиваете, если ваш код DCT является правильным или что-то?
Джим Клэй
Спасибо за ваши Коментарии. Я добавил еще одно предложение в начале. Моя цель - реализовать ПКТ на основе БПФ.
Framester

Ответы:

18

Я читал об этом , и есть несколько способов сделать это, используя разного размера N. Мой Matlab ржавый, так что здесь они находятся в Python ( Nэто длина входного сигнала x, kявляется arange(N)= ):[0,1,2,...,N1]

DCT типа 2 с использованием 4N FFT и без смен

Сигнал [a, b, c, d]становится

[0, a, 0, b, 0, c, 0, d, 0, d, 0, c, 0, b, 0, a],

Затем возьмите БПФ, чтобы получить спектр

[A, B, C, D, 0, -D, -C, -B, -A, -B, -C, -D, 0, D, C, B]

затем выбросьте все, кроме первого [A, B, C, D], и все готово

u = zeros(4 * N)
u[1:2*N:2] = x
u[2*N+1::2] = x[::-1]

U = fft(u)[:N]
return U.real

Тип 2 DCT с использованием 2N БПФ зеркальный (Махул)

[a, b, c, d]становится [a, b, c, d, d, c, b, a]. Возьмите БПФ этого, чтобы получить [A, B, C, D, 0, D*, C*, B*], затем выбросьте все, но [A, B, C, D]и умножьте это на (сдвиг половины выборкиejπk2N ), чтобы получить DCT:

y = empty(2*N)
y[:N] = x
y[N:] = x[::-1]

Y = fft(y)[:N]

Y *= exp(-1j*pi*k/(2*N))
return Y.real

Тип 2 DCT с использованием 2N FFT (Makhoul)

[a, b, c, d]становится [a, b, c, d, 0, 0, 0, 0]. Возьмите БПФ этого, чтобы получить [A, B, C, D, E, D*, C*, B*], затем выбросьте все, но [A, B, C, D]и умножьте это на чтобы получить DCT:2ejπk2N

y = zeros(2*N)
y[:N] = x

Y = fft(y)[:N]

Y *= 2 * exp(-1j*pi*k/(2*N))
return Y.real

DCT типа 2 с использованием N FFT (Makhoul)

Сигнал [a, b, c, d, e, f]становится [a, c, e, f, d, b], затем возьмите БПФ, чтобы получить [A, B, C, D, C*, B*]. Умножьте на а затем принять реальное участие:2ejπk2N

v = empty_like(x)
v[:(N-1)//2+1] = x[::2]

if N % 2: # odd length
    v[(N-1)//2+1:] = x[-2::-2]
else: # even length
    v[(N-1)//2+1:] = x[::-2]

V = fft(v)

V *= 2 * exp(-1j*pi*k/(2*N))
return V.real

На моей машине они все примерно одинаковой скорости, поскольку генерация exp(-1j*pi*k/(2*N))занимает больше времени, чем FFT. : D

In [99]: timeit dct2_4nfft(a)
10 loops, best of 3: 23.6 ms per loop

In [100]: timeit dct2_2nfft_1(a)
10 loops, best of 3: 20.1 ms per loop

In [101]: timeit dct2_2nfft_2(a)
10 loops, best of 3: 20.8 ms per loop

In [102]: timeit dct2_nfft(a)
100 loops, best of 3: 16.4 ms per loop

In [103]: timeit scipy.fftpack.dct(a, 2)
100 loops, best of 3: 3 ms per loop
эндолиты
источник
2
Отличный ответ, очень помог с моей реализацией! Дополнительное примечание: последний метод «DCT типа 2 с использованием N FFT» по-прежнему работает должным образом, если длина сигнала нечетна; последний элемент перемещается к среднему элементу. Я проверил математику и код для этого факта.
Наюки,
1
@ Наюки Вы генерируете exp(-1j*pi*k/(2*N))или есть ярлык для этого шага?
эндолиты
Я генерирую exp(-1j*pi*k/(2*N))в своем коде , потому что сдвиг четверти выборки необходим для работы отображения DCT-DFT. Что заставляет вас спрашивать?
Наюки
Привет, как это будет работать для DCT типа III, чтобы вычислить обратное DCT-II?
Джек Х
8

Я не уверен, какой метод использует бумага, на которую вы ссылались, но я вывел это раньше, и вот что я получил в итоге: (Предполагая, что - это ввод)x(n)

позволять

y(n)={x(n),n=0,1,...,N1x(2N1n),n=N,N+1,...,2N1

DCT тогда дается

C(k)=Re{ejπk2NFFT{y(n)}}

Таким образом, вы в основном создаете последовательность длиной y ( n ), где первая половина - x ( n ), а вторая половина - x ( n ).2Ny(n)x(n)x(n) обратном порядке. Затем просто возьмите БПФ и умножьте этот вектор на фазовый сдвиг. Наконец, возьмите только реальную часть, и у вас есть DCT.

Вот код в MATLAB.

function C = fdct(x)
    N = length(x);
    y = zeros(1,2*N);
    y(1:N) = x;
    y(N+1:2*N) = fliplr(x);
    Y = fft(y);
    k=0:N-1;
    C = real(exp(-j.* pi.*k./(2*N)).*Y(1:N));

Редактировать:

Примечание. Используемая формула DCT:

C(k)=2n=0N1x(n)cos(πk2N(2n+1))

Существует несколько способов масштабирования суммирования, поэтому оно может не совпадать точно с другими реализациями. Например, MATLAB использует:

C(k)=w(k)n=0N1x(n)cos(πk2N(2n+1))

w(0)=1Nw(1...N1)=2N

Вы можете объяснить это, правильно масштабируя вывод.

Джейсон Б
источник
1
у (n) должна быть N-длины, а не 2N-длины. Таким образом, вы получаете 4-кратную скорость вычислений, вычисляя DCT N-длины из сигнала N-длины вместо 2N FFT из 2N-сигнала. fourier.eng.hmc.edu/e161/lectures/dct/node2.html www-ee.uta.edu/dip/Courses/EE5355/Discrete%20class%201.pdf
эндолиты
0

Для истинных научных вычислений также важно количество используемой памяти. Поэтому точка N БПФ более привлекательна для меня. Это возможно только благодаря эрмитовой симметрии сигнала. Ссылка Махул дана здесь. И фактически имеет алгоритм расчета DCT и IDCT или DCT10 и DCT01.
http://ieeexplore.ieee.org/abstract/document/1163351/

Hasbestein
источник