Дискретная косинусная трансформация 1D Matlab

Question

Я пытаюсь реализовать DCT в MATLAB с помощью умножения матрицы-вектора. В частности, я хотел бы создать DCT-матрицу коэффициентов, а затем использовать ее для умножения на 1D-сигнал для вычисления 1D DCT.

Вот мой код, чтобы получить матрицу DCT:

function D=dct1d(n) 
for u=0:n-1 
    if u==0 
     au=sqrt(1/n); 
    else 
     au=sqrt(2/n); 
    end 
    for x=0:n-1 
     D(u+1,x+1)=au*cos(((2*x+1)*u*pi)/2*n); 
    end 
end 

После этого я попытался выполнения DCT с тестовым сигналом x = [1 2 3 4 5 6 7 8]:

x=[1 2 3 4 5 6 7 8]; 
y=dct1(8)*x'; 

Ответ дает это:

12.7279 
18.0000 
18.0000 
18.0000 
18.0000 
18.0000 
18.0000 
18.0000 

Однако правильный ответ:

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

Answer 1

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

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); 
%//        ^^^ 

Because multiplication and division are equal in precedence, это точно так же, как делают:

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. BTW, один из способов проверить, есть ли у вас правильный ответ, использовать функцию dctmtx, которая вычисляет матрицу коэффициентов DCT 10, которую вы ищете. Однако эта функция является частью панели обработки сигналов, поэтому, если у вас ее нет, вы, к сожалению, не сможете использовать эту функцию, но если я могу предложить альтернативный ответ вместо использования циклов 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 под капотом, после того как вы удалите все проверки ошибок и проверки согласованности ввода.}