Transformation de cosinus discrète 1D Matlab

Question

J'essaie de mettre en œuvre DCT dans MATLAB en utilisant la multiplication matrice-vecteur. Spécifiquement, je voudrais créer la matrice de coefficients DCT puis l'utiliser pour multiplier avec le signal 1D pour calculer le DCT 1D.

Voici mon code pour produire la matrice TCD:

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 

Après cela, j'ai essayé d'effectuer le TCD avec un signal de test de x = [1 2 3 4 5 6 7 8]:

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

La réponse qu'il donne est:

12.7279 
18.0000 
18.0000 
18.0000 
18.0000 
18.0000 
18.0000 
18.0000 

Cependant, la bonne réponse est:

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

Answer 1

L'erreur dans votre code est très légère mais fondamentale. La ligne où vous calculer les coefficients:

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

Jetez un oeil à la dernière partie de la ligne:

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

Because multiplication and division are equal in precedence, c'est exactement la même chose que faire:

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

En tant que tel, vous ne divisez pas par 2n. Vous divisez par 2 puis en multipliant par n ce qui est incorrect. Il vous suffit d'entourer l'opération 2*n avec des supports:

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

Une fois que vous faites cela, nous obtenons la matrice correcte TCD. BTW, une façon de vérifier si vous avez la bonne réponse est d'utiliser la fonction dctmtx, qui calcule la matrice de coefficients DCT N x N que vous recherchez. Cependant, cette fonction fait partie de Signal Processing Toolbox, donc si vous ne l'avez pas, alors vous ne pouvez malheureusement pas utiliser cette fonction, mais si je peux suggérer une réponse alternative au lieu d'utiliser les boucles for, je créerais une grille 2D des coordonnées avec meshgrid, puis calculez les produits par éléments.

Quelque chose comme ça fonctionnerait à la place:

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 

Au lieu d'utiliser une déclaration if pour déterminer ce poids, nous devons appliquer par ligne, nous pouvons simplement utiliser sqrt(2/n) puis diviser par sqrt(2) pour la première ligne de telle sorte que vous divisez par sqrt(1/n) à la place. Ce code devrait produire les mêmes résultats que votre code corrigé .

---

Dans tous les cas, une fois que je fait ces corrections, j'ai comparé les deux réponses entre ce que votre code donne et ce que dctmtx donne et il a raison:

>> 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 

Une fois que nous obtenons la matrice TCD corrigée , nous pouvons vérifier les calculs TCD réels en effectuant la multiplication de la matrice avec le vecteur de test de 1:8 que vous avez utilisé:

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

ans = 

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

---

^{1. Ce code est en fait ce qui est fait dans dctmtx sous le capot, une fois que vous supprimez tous les contrôles d'erreurs et de cohérence d'entrée.}