Дискретное косинусное преобразование 1D Matlab

Ошибка в вашем коде очень небольшая, но фундаментальная. Строка, в которой вы вычисляете коэффициенты:

D(u+1,x+1)=au*cos(((2*x+1)*u*pi)/2*n);

Взгляните на самую последнюю часть строки:

D(u+1,x+1)=au*cos(((2*x+1)*u*pi)/2*n);
%//                              ^^^

Поскольку умножение и деление имеют одинаковый приоритет, это именно тот то же, что делать:

D(u+1,x+1)=au*cos((((2*x+1)*u*pi)/2)*n);

Таким образом, вы не делите на 2n. Вы делите на 2, а затем умножаете на n, что неверно. Вам просто нужно заключить операцию 2*n в скобки:

D(u+1,x+1)=au*cos(((2*x+1)*u*pi)/(2*n)); 

Как только вы это сделаете, мы получим правильную матрицу DCT. Кстати, один из способов проверить правильность ответа — использовать dctmtx, которая вычисляет искомую матрицу коэффициентов N x N DCT. Тем не менее, эта функция является частью набора инструментов обработки сигналов, поэтому, если у вас ее нет, вы, к сожалению, не можете использовать эту функцию, но если я могу предложить альтернативный ответ вместо использования for циклов, я бы построил 2D-сетку. координат с помощью meshgrid, а затем вычислить поэлементные произведения.

Вместо этого будет работать что-то вроде этого:

function D = dct1d(n)

[x,u] = meshgrid(0:n-1);
D = sqrt(2/n)*cos(((2*x+1).*u*pi)/(2*n)); 
D(1,:) = D(1,:) / sqrt(2);

end

Вместо использования оператора if для определения веса, который нам нужно применить к строке, мы можем просто использовать sqrt(2/n), а затем разделить на sqrt(2) для первой строки, чтобы вместо этого вы делили на sqrt(1/n). Этот код должен дать те же результаты, что и исправленный код¹.

В любом случае, как только я внес эти исправления, я сравнил оба ответа между тем, что дает ваш код, и тем, что дает dctmtx, и это правильно:

>> dct1d(8)

ans =

    0.3536    0.3536    0.3536    0.3536    0.3536    0.3536    0.3536    0.3536
    0.4904    0.4157    0.2778    0.0975   -0.0975   -0.2778   -0.4157   -0.4904
    0.4619    0.1913   -0.1913   -0.4619   -0.4619   -0.1913    0.1913    0.4619
    0.4157   -0.0975   -0.4904   -0.2778    0.2778    0.4904    0.0975   -0.4157
    0.3536   -0.3536   -0.3536    0.3536    0.3536   -0.3536   -0.3536    0.3536
    0.2778   -0.4904    0.0975    0.4157   -0.4157   -0.0975    0.4904   -0.2778
    0.1913   -0.4619    0.4619   -0.1913   -0.1913    0.4619   -0.4619    0.1913
    0.0975   -0.2778    0.4157   -0.4904    0.4904   -0.4157    0.2778   -0.0975

>> dctmtx(8)

ans =

    0.3536    0.3536    0.3536    0.3536    0.3536    0.3536    0.3536    0.3536
    0.4904    0.4157    0.2778    0.0975   -0.0975   -0.2778   -0.4157   -0.4904
    0.4619    0.1913   -0.1913   -0.4619   -0.4619   -0.1913    0.1913    0.4619
    0.4157   -0.0975   -0.4904   -0.2778    0.2778    0.4904    0.0975   -0.4157
    0.3536   -0.3536   -0.3536    0.3536    0.3536   -0.3536   -0.3536    0.3536
    0.2778   -0.4904    0.0975    0.4157   -0.4157   -0.0975    0.4904   -0.2778
    0.1913   -0.4619    0.4619   -0.1913   -0.1913    0.4619   -0.4619    0.1913
    0.0975   -0.2778    0.4157   -0.4904    0.4904   -0.4157    0.2778   -0.0975

Как только мы получим скорректированную матрицу DCT, мы можем проверить фактические вычисления DCT, выполнив умножение матрицы на тестовый вектор 1:8, который вы использовали:

>> dct1d(8)*((1:8).')

ans =

   12.7279
   -6.4423
   -0.0000
   -0.6735
         0
   -0.2009
   -0.0000
   -0.0507

^{1. Этот код на самом деле делается в dctmtx под капотом, как только вы уберете все проверки ошибок и проверки согласованности ввода.}