Я ищу простой алгоритм для быстрого DCT и IDCT матрицы [NxM]

Question

Я ищу простой алгоритм для выполнения быстрого DCT (тип 2) матрицы любого размера [NxM], а также алгоритм для обратного преобразования IDCT (также называемый DCT тип 3).

Мне нужен алгоритм DCT-2D, но даже алгоритм DCT-1D достаточно хорош, потому что я могу использовать DCT-1D для реализации DCT-2D (и IDCT-1D для реализации IDCT-2D).

PHP-код предпочтительнее, но подойдет любой достаточно понятный алгоритм.

Мой текущий PHP-скрипт для реализации DCT/IDCT очень медленный, когда размер матрицы превышает [200x200].

Я прыгал, чтобы найти способ преформировать DCT до [4000x4000] менее чем за 20 секунд. кто нибудь знает как это сделать?

Answer 1

Вот мое вычисление 1D FDCT и IFDCT для FFT с одинаковой длиной:

//---------------------------------------------------------------------------
void DFCTrr(double *dst,double *src,double *tmp,int n)
    {
    // exact normalized DCT II by N DFFT
    int i,j;
    double nn=n,a,da=(M_PI*(nn-0.5))/nn,a0,b0,a1,b1,m;
    for (j=  0,i=n-1;i>=0;i-=2,j++) dst[j]=src[i];
    for (j=n-1,i=n-2;i>=0;i-=2,j--) dst[j]=src[i];
    DFFTcr(tmp,dst,n);
    m=2.0*sqrt(2.0);
    for (a=0.0,j=0,i=0;i<n;i++,j+=2,a+=da)
        {
        a0=tmp[j+0]; a1= cos(a);
        b0=tmp[j+1]; b1=-sin(a);
        a0=(a0*a1)-(b0*b1);
        if (i) a0*=m; else a0*=2.0;
        dst[i]=a0;
        }
    }
//---------------------------------------------------------------------------
void iDFCTrr(double *dst,double *src,double *tmp,int n)
    {
    // exact normalized DCT III = iDCT II by N iDFFT
    int i,j;
    double nn=n,a,da=(M_PI*(nn-0.5))/nn,a0,m,aa,bb;
    m=1.0/sqrt(2.0);
    for (a=0.0,j=0,i=0;i<n;i++,j+=2,a+=da)
        {
        a0=src[i];
        if (i) a0*=m;
        aa= cos(a)*a0;
        bb=+sin(a)*a0;
        tmp[j+0]=aa;
        tmp[j+1]=bb;
        }
    m=src[0]*0.25;
    iDFFTrc(src,tmp,n);
    for (j=  0,i=n-1;i>=0;i-=2,j++) dst[i]=src[j]-m;
    for (j=n-1,i=n-2;i>=0;i-=2,j--) dst[i]=src[j]-m;
    }
//---------------------------------------------------------------------------

• dst вектор назначения [n]
• src является исходным вектором [n]
• tmp является временным вектором [2n]

Эти массивы не должны перекрываться!!! Это взято из моего класса преобразования, так что я надеюсь, что не забыл что-то скопировать.

• XXXrr означает, что назначение является реальным, а источник - также реальным доменом
• XXXrc означает, что назначение является реальным, а источник - сложным
• XXXcr означает, что назначение является сложным, а источник - реальным доменом

Все данные double для массивов, для комплексного домена первое число является действительным, а второе - мнимой частью, поэтому массив 2N размер. Обе функции используют FFT и iFFT, если вам нужен код, также для них прокомментируйте меня. Просто чтобы быть уверенным, что я добавил не быстрое внедрение их ниже. Это гораздо проще скопировать, потому что быстрые используют слишком много иерархии классов преобразования.

медленные реализации DFT, iDFT для тестирования:

//---------------------------------------------------------------------------
void transform::DFTcr(double *dst,double *src,int n)
    {
    int i,j;
    double a,b,a0,_n,q,qq,dq;
    dq=+2.0*M_PI/double(n); _n=2.0/double(n);
    for (q=0.0,j=0;j<n;j++,q+=dq)
        {
        a=0.0; b=0.0;
        for (qq=0.0,i=0;i<n;i++,qq+=q)
            {
            a0=src[i];
            a+=a0*cos(qq);
            b+=a0*sin(qq);
            }
        dst[j+j  ]=a*_n;
        dst[j+j+1]=b*_n;
        }
    }
//---------------------------------------------------------------------------
void transform::iDFTrc(double *dst,double *src,int n)
    {
    int i,j;
    double a,a0,a1,b0,b1,q,qq,dq;
    dq=+2.0*M_PI/double(n);
    for (q=0.0,j=0;j<n;j++,q+=dq)
        {
        a=0.0;
        for (qq=0.0,i=0;i<n;i++,qq+=q)
            {
            a0=src[i+i  ]; a1=+cos(qq);
            b0=src[i+i+1]; b1=-sin(qq);
            a+=(a0*a1)-(b0*b1);
            }
        dst[j]=a*0.5;
        }
    }
//---------------------------------------------------------------------------

Так что для тестирования просто переписать имена DFFTcr а также iDFFTrc (или используйте их для сравнения с вашим FFT,iFFT) когда код работает должным образом, реализуйте свой собственный FFT, iFFT. Подробнее об этом см.:

• Как вычислить дискретное преобразование Фурье?

2D DFCT

1. изменить размер src матрица к власти 2

   добавляя нули, чтобы использовать быстрый алгоритм, размер должен всегда быть степенью 2!!!

2. выделять NxN реальные матрицы tmp,dst а также 1xN сложный вектор t

3. преобразовать строки DFCTrr

   DFCT(tmp.line(i),src.line(i),t,N)
   

4. транспонировать tmp матрица

5. преобразовать строки DFCTrr

   DFCT(dst.line(i),tmp.line(i),t,N)
   

6. транспонировать dst матрица

7. нормализовать dst умножить матрицу на 0.0625

2D iDFCT

То же, что и выше, но использовать iDFCTrr и умножить на 16.0 вместо.

[Заметки]

Перед реализацией ваших собственных FFT и iFFT убедитесь, что они дают тот же результат, что и мой, иначе DCT/iDCT не будет работать должным образом!!!