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

Я пытаюсь реализовать 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
1 3
matlab image-compression signal-processing

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

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

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) вместо этого. Этот код должен давать те же результаты, что и ваш исправленный код1.

В любом случае, как только я сделал эти исправления, я сравнил оба ответа между тем, что дает ваш код, и тем, что дает 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 под капотом, как только вы убираете все проверки ошибок и проверки согласованности ввода.

3