DCT C语言 修改了c程序以及matlab的程序

DCT

定义

给定序列x(n),n=0,1…N-1,其离散余弦变换定义为

$X_c(0)=\frac{1}{\sqrt{N}}\sum_{n=0}^{N-1}x(n)$Xc​(0)=N​1​∑n=0N−1​x(n)

$X_c(k)=\sqrt{\frac{2}{N}}\sum_{n=0}^{N-1}x(n)cos\frac{(2n+1)*k*\pi}{2N},n=0,1...,N-1$Xc​(k)=N2​​∑n=0N−1​x(n)cos2N(2n+1)∗k∗π​,n=0,1...,N−1

其变换的核函数

$C_{k,n}=\sqrt{\frac{2}{N}}g_kcos\frac{(2n+1)k\pi}{2N}, k,n=0,1....,N-1$Ck,n​=N2​​gk​cos2N(2n+1)kπ​,k,n=0,1....,N−1

是实数，式中系数

c程序

#include<stdio.h>
#include<math.h>
#define N 64
#define pi 3.14159
int main()
{
    int n,k;
    float y[100],xn[100];//输出，输入信号
    for(k=0; k<N; k++)//循环求X(k)
    {
        if(k==0)
        {
            for(n=0; n<N; n++)
            {
                xn[n]=cos(n*pi/6.0);
                y[n]=sqrt(1.0/N)*xn[n];//公式
                y[k]+=y[n];
            }

        }
        else
        {
            for(n=0; n<N; n++)
            {
                xn[n]=cos(n*pi/6.0);
                y[n]=sqrt(2.0/N)*xn[n]*cos((2.0*n+1)*k*pi/(2.0*N));//公式
                y[k]+=y[n];
            }
        }
        printf(" %d   %lf\n",k,y[k]);
    }

}

matlab

>> N=64;
>> n=0:1:N;
>> xn=cos(n*pi/6);
>> figure
>> y=dct(xn);
>> stem(n,y)

c结果

matlab结果

对比有点差异，跟FFT的差异相似。
有大佬知道怎么修改的请评论下，感谢

测试的问题

#include <stdio.h>//标准输入输出头文
#include<math.h>//声明了常用的一些数学运算
#define pi 3.1415
#define N 64//64个点
int main()
{
    int n,k;
    double xn;//输入信号
    double y[100];//输出，输入信号
    for(k=0; k<N; k++)//循环求X(k)
    {
        for(n=0; n<N; n++)
        {

                xn=1.0;
                y[k] += xn;

        }
        printf("%d  %f\n",k,y[k]);
    }

}

得出的数据

0
-1.#QNAN0
1 64.000000
2
-1.#QNAN0
3 64.000000
4 64.000000
5 64.000000
6 64.000000
7
-1.#QNAN0
8 64.000000
9 64.000000
10 3072570228342306800000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000
11 64.000000
12 64.000000
13 64.000000
14 64.000000
15 3072569644077048600000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000
16 64.000000
17 64.000000
18 64.000000
19 64.000000
20 64.000000
21 64.000000
22 64.000000
23 64.000000
24 64.000000
25 64.000000
26 64.000000
27 64.000000
28 64.000000
29 64.000000
30 64.000000
31 2833760479163969400000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000
32 64.000000
33 64.000000
34 64.000000
35 64.000000
36 64.000000
37 64.000000
38 64.000000
39 64.000000
40 64.000000
41 3488645024638201700000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000
42
-1.#QNAN0
43 3581851895918701500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000
44 64.000000
45 64.000000
46 64.000000
47 64.000000
48 64.000000
49 64.000000
50 64.000000
51 64.000000
52 64.000000
53 64.000000
54 64.000000
55 64.000000
56 64.000000
57 64.000000
58 64.000000
59 64.000000
60 64.000000
61 64.000000
62 64.000000
63 64.000000

发现 0，2，7，10，15，31，41，42，43是错误的。

修改成功后

发现是定义的范围大小不够,以及增添一个下x[n]函数，改了如下

double y[1000]，x[100]

#include<stdio.h>
#include<math.h>
#define N 64
#define pi 3.14159
int main()
{
    int n,k;
    double xn;
    double  y[1000],x[100];//输出，输入信号
    for(k=0; k<N; k++)//循环求X(k)
    {
        if(k==0)
        {
            for(n=0; n<N; n++)
            {
                xn=cos(n*pi/6.0);
                x[n]=sqrt(1.0/N)*xn;//公式
                y[k]+=x[n];
            }

        }
        else
        {
            for(n=0; n<N; n++)
            {
                xn=cos(n*pi/6.0);
                x[n]=sqrt(2.0/N)*xn*cos((2.0*n+1)*k*pi/(2.0*N));//公式
                y[k]+=x[n];
            }
        }
        printf(" %d   %lf\n",k,y[k]);
    }

}

然后我发现matlab的好像程序错了，修改了下

> n=0:1:63;
>> x=cos(n*pi/6);
>> y=dct(x);
>> figure
>> stem(n,y)

c 的结果图

matlab的结果图

结果图跟matlab相差不多