C语言实现离散余弦变换(DCT)并用MATLAB和Python验证
概念
離散余弦變換(Discrete Cosine Transform,DCT)是可分離的變換,其變換核為余弦函數。是與傅里葉變換相關的一種變換,它相當于把離散傅里葉變換的虛數部分丟掉,只使用實數。DCT除了具有一般的正交變換性質外, 它的變換陣的基向量能很好地描述人類語音信號和圖像信號的相關特征。因此,在對語音信號、圖像信號的變換中,DCT變換被認為是一種準最佳變換。
設x(n)是N個有限值的一維實數信號序列,n=0,1,?,N?1n=0,1,\cdots,N-1n=0,1,?,N?1的完備正交歸一函數如下:
X(k)=a(k)∑n=0N?1x(n)cos?[π(2n+1)k2N]x(n)=a(k)∑k=0N?1X(k)cos?[π(2n+1)k2N]\begin{aligned} X(k) &=a(k) \sum_{n=0}^{N-1} x(n) \cos \left[\frac{\pi(2 n+1) k}{2 N}\right] \\ x(n) &=a(k) \sum_{k=0}^{N-1} X(k) \cos \left[\frac{\pi(2 n+1) k}{2 N}\right] \end{aligned} X(k)x(n)?=a(k)n=0∑N?1?x(n)cos[2Nπ(2n+1)k?]=a(k)k=0∑N?1?X(k)cos[2Nπ(2n+1)k?]?
其中a(k)a(k)a(k)定義為:
a(k)={1/Nk=02/N1≤k≤N?1a(k)=\left\{\begin{array}{ll} \sqrt{1 / N} & k=0 \\ \sqrt{2 / N} & 1 \leq k \leq N-1 \end{array}\right. a(k)={1/N?2/N??k=01≤k≤N?1?
由上述兩個式子可以得到DCT的另一種表現形式:
X(k)=2N∑n=0NC(k)x(n)cos?[π(2π+1)k2N]k=0,1…,N?1X(k)=\sqrt{\frac{2}{N}} \sum_{n=0}^{N} C(k) x(n) \cos \left[\frac{\pi(2 \pi+1) k}{2 N}\right] \quad k=0,1 \ldots, N-1 X(k)=N2??n=0∑N?C(k)x(n)cos[2Nπ(2π+1)k?]k=0,1…,N?1
任務需求
我們產生一個幅度為1,頻率為100hz,采樣頻率為5000的正弦波信號,并對其進行DFC變換,將產生的正弦波信號保存在text.txt中,將經過DCT變換的數據保存在result.txt文件中。再把text.txt數據放到MATLAB和Python中驗證,檢驗數據計算是否正確。
C語言實現
在C語言環境下的代碼如下:
#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就是長度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++){//產生信號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) //計算k=0時的akak=sqrt(1.0/len); else //計算k!=0時的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]等于累和結果s乘以系數C}for(k=0;k<len;k++) {//printf("%dDFC:%lf\n",k,X[k]);fprintf(fp2, "%.4f\n", X[k]);//輸出離散余弦變換結果 }printf("successfully!"); }在cmd中輸入如下指令進行編譯和運行
gcc dct-c.c a.exe 1 100 5000
打開test.txt和result.txt文件,查看輸出的數據是否是5000行,輸出結果如下。
在gnuplot中輸入如下的代碼繪制出圖像
MATLAB驗證
在MATLAB環境下的代碼如下:
close all; x=load("test.txt");%加載數據 subplot(211);%兩個子圖 plot(x)%繪制原始的信號圖 xlim([0,500]);%控制x軸范圍 y=dct(x);%進行dct變換 subplot(212);%第二個子圖 plot(y);%繪制變換后的頻譜圖MATLAB中的輸出結果
Python實現
在Python環境下的代碼如下:
import numpy as np from scipy.fftpack import dct import matplotlib.pyplot as plt file = 'test.txt' data = np.loadtxt(file)#讀入數據 a = len(data)#計算數據長度 x1 = np.linspace(0, 1, a)#創建0-1的等差數組 a個 y = dct(data)#對數據進行dct變換 plt.subplot(211)#繪制第一幅圖 plt.xlim((0, 0.1))#控制xy軸范圍 plt.ylim((-2, 2))#圖表標題 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中的輸出結果如圖
可以驗證我們C語言實現的算法是正確的。
總結
以上是生活随笔為你收集整理的C语言实现离散余弦变换(DCT)并用MATLAB和Python验证的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 怎么从mysql注册表删除用户_mysq
- 下一篇: java 上下文加载器_【深入理解Jav