概念

離散余弦變換（Discrete Cosine Transform，DCT）是可分離的變換，其變換核為余弦函數(shù)。是與傅里葉變換相關(guān)的一種變換，它相當(dāng)于把離散傅里葉變換的虛數(shù)部分丟掉，只使用實(shí)數(shù)。DCT除了具有一般的正交變換性質(zhì)外， 它的變換陣的基向量能很好地描述人類語(yǔ)音信號(hào)和圖像信號(hào)的相關(guān)特征。因此，在對(duì)語(yǔ)音信號(hào)、圖像信號(hào)的變換中，DCT變換被認(rèn)為是一種準(zhǔn)最佳變換。
設(shè)x(n)是N個(gè)有限值的一維實(shí)數(shù)信號(hào)序列，$n=0,1,?,N?1$的完備正交歸一函數(shù)如下：
$$\begin{array}{rl}X(k)& =a(k)\sum _{n=0}^{N?1}x(n)\cos \left[\frac{\pi (2n+1)k}{2N}\right]\\ x(n)& =a(k)\sum _{k=0}^{N?1}X(k)\cos \left[\frac{\pi (2n+1)k}{2N}\right]\end{array}$$
其中$a(k)$定義為：
$$a(k)=\{\begin{array}{ll}\sqrt{1/N}& k=0\\ \sqrt{2/N}& 1\le k\le N?1\end{array}$$
由上述兩個(gè)式子可以得到DCT的另一種表現(xiàn)形式：
$$X(k)=\sqrt{\frac{2}{N}}\sum _{n=0}^{N}C(k)x(n)\cos \left[\frac{\pi (2\pi +1)k}{2N}\right]k=0,1\dots ,N?1$$

任務(wù)需求

我們產(chǎn)生一個(gè)幅度為1，頻率為100hz，采樣頻率為5000的正弦波信號(hào)，并對(duì)其進(jìn)行DFC變換，將產(chǎn)生的正弦波信號(hào)保存在text.txt中，將經(jīng)過(guò)DCT變換的數(shù)據(jù)保存在result.txt文件中。再把text.txt數(shù)據(jù)放到MATLAB和Python中驗(yàn)證,檢驗(yàn)數(shù)據(jù)計(jì)算是否正確。

C語(yǔ)言實(shí)現(xiàn)

在C語(yǔ)言環(huán)境下的代碼如下：

#include <stdio.h> #include<math.h> #include<stdlib.h> #define PI 3.1415926 void main(int argc,char* argv[]) {double A = atof(argv[1]);//幅度double f = atof(argv[2]);//頻率/hzdouble fn = atof(argv[3]);//采樣頻率int len = (int)fn;//ken就是長(zhǎng)度Ndouble x[len];double X[len];double t,data;int i,j,k;double ak;FILE* fp1 = fopen("test.txt","w");FILE* fp2 = fopen("result.txt","w");for (i = 0; i < len; i++){//產(chǎn)生信號(hào)t = i / fn;x[i] = A * sin(2 * PI * f * t);//printf("%d %lf\n", i, x[i]);fprintf(fp1, "%lf\n", x[i]);}//DFT變換//求出a(k)for(k=0;k<len;k++){double data =0;if(k==0) //計(jì)算k=0時(shí)的akak=sqrt(1.0/len); else //計(jì)算k!=0時(shí)的akak=sqrt(2.0/len);for(j=0;j<len;j++)//離散余弦變換 {t = x[j]*cos(((2*j+1)*k*PI)/(2*len));data = data + t; //累加求和}X[k]=ak*data;//X[k]等于累和結(jié)果s乘以系數(shù)C}for(k=0;k<len;k++) {//printf("%dDFC:%lf\n",k,X[k]);fprintf(fp2, "%.4f\n", X[k]);//輸出離散余弦變換結(jié)果 }printf("successfully!"); }

在cmd中輸入如下指令進(jìn)行編譯和運(yùn)行

gcc dct-c.c a.exe 1 100 5000

打開test.txt和result.txt文件，查看輸出的數(shù)據(jù)是否是5000行，輸出結(jié)果如下。

在gnuplot中輸入如下的代碼繪制出圖像

set multiplot layout 2,1 set xrange[0:500] plot "test.txt" w l unset xrange plot "result.txt" w l

MATLAB驗(yàn)證

在MATLAB環(huán)境下的代碼如下：

close all; x=load("test.txt");%加載數(shù)據(jù) subplot(211);%兩個(gè)子圖 plot(x)%繪制原始的信號(hào)圖 xlim([0,500]);%控制x軸范圍 y=dct(x);%進(jìn)行dct變換 subplot(212);%第二個(gè)子圖 plot(y);%繪制變換后的頻譜圖

MATLAB中的輸出結(jié)果

Python實(shí)現(xiàn)

在Python環(huán)境下的代碼如下：

import numpy as np from scipy.fftpack import dct import matplotlib.pyplot as plt file = 'test.txt' data = np.loadtxt(file)#讀入數(shù)據(jù) a = len(data)#計(jì)算數(shù)據(jù)長(zhǎng)度 x1 = np.linspace(0, 1, a)#創(chuàng)建0-1的等差數(shù)組 a個(gè) y = dct(data)#對(duì)數(shù)據(jù)進(jìn)行dct變換 plt.subplot(211)#繪制第一幅圖 plt.xlim((0, 0.1))#控制xy軸范圍 plt.ylim((-2, 2))#圖表標(biāo)題 plt.title('signal', fontsize=12, color='black') plt.plot(x1,data)plt.subplot(212)#繪制第二幅圖 plt.title('dct', fontsize=12, color='black') plt.plot(y)#控制兩幅子圖的上下間距 plt.subplots_adjust(left=None, bottom=None, right=None, top=None,wspace=None, hspace=1) plt.savefig('1.7.eps')#保存為eps矢量圖 plt.show()

在Python中的輸出結(jié)果如圖
可以驗(yàn)證我們C語(yǔ)言實(shí)現(xiàn)的算法是正確的。

總結(jié)

以上是生活随笔為你收集整理的C语言实现离散余弦变换(DCT)并用MATLAB和Python验证的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。