-
离散余弦变换原理
一句话概括核心原理：频率不同的正弦波相乘，其黎曼积分总是为零。下图是网上找来的离散余弦变换公式：
-
C代码实现
dct.h
//==============================================================================
// Copyright (C) 2019 王小康. All rights reserved.
//
// 作者: 王小康
// 描述: 一维二维DCT变换
// 日期: 2019-05-09
//
//==============================================================================
#ifndef _DCT_H_
#define _DCT_H_
#define DCT_LEN 8 //一维或二维矩阵边长
void dct_init (void);
void dct_1Ddct (float *out, float *in);
void dct_1Didct (float *out, float *in);
void dct_1Dspectrum (float *out, float *in);
void dct_2Ddct (float out[][DCT_LEN], float in[][DCT_LEN]);
void dct_2Didct (float out[][DCT_LEN], float in[][DCT_LEN]);
#endif
dct.c
//==============================================================================
// Copyright (C) 2019 王小康. All rights reserved.
//
// 作者: 王小康
// 描述: 一维二维DCT变换
// 日期: 2019-05-09
//
//==============================================================================
#include <math.h>
#include "dct.h"
static float dctTableA[DCT_LEN][DCT_LEN];
static float sqrt_1_N;
static float sqrt_2_N;
////////////////////////////////////////////////////////////////////////////////
// 离散余弦变换初始化，生成变换核。一维二维同核。
void dct_init(void){
sqrt_1_N = sqrt(1.0/DCT_LEN);
sqrt_2_N = sqrt(2.0/DCT_LEN);
int i, j;
for(i=0; i<DCT_LEN; i++){
dctTableA[0][i] = sqrt_1_N;
}
for(i=1; i<DCT_LEN; i++){
for(j=0; j<DCT_LEN; j++){
dctTableA[i][j] = sqrt_2_N * cos((3.14159265359*i*(2*j+1))/(2*DCT_LEN));
}
}
}
// 一维离散余弦变换
void dct_1Ddct(float *out, float *in){
unsigned i, j;
for(i=0; i<DCT_LEN; i++){
out[i] = 0;
for(j=0; j<DCT_LEN; j++) out[i] += dctTableA[i][j] * in[j];
}
}
// 一维离散余弦变换，返回频谱图
void dct_1Dspectrum(float *out, float *in){
unsigned i, j;
out[0] = 0;
for(i=0; i<DCT_LEN; i++) out[0] += in[i];
out[0] /= DCT_LEN;
if(out[0] < 0) out[0] = -out[0];
for(i=1; i<DCT_LEN; i++){
out[i] = 0;
for(j=0; j<DCT_LEN; j++) out[i] += dctTableA[i][j] * in[j];
out[i] *= sqrt_2_N;
if(out[i] < 0) out[i] = -out[i];
}
}
// 一维离散余弦逆变换
void dct_1Didct(float *out, float *in){
unsigned i, j;
for(i=0; i<DCT_LEN; i++){
out[i] = 0;
for(j=0; j<DCT_LEN; j++) out[i] += dctTableA[j][i] * in[j];
}
}
// 二维离散余弦变换
void dct_2Ddct(float out[][DCT_LEN], float in[][DCT_LEN]){
float tmp[DCT_LEN][DCT_LEN];
float sum;
int i, j, k;
for(i=0; i<DCT_LEN; i++){
for(j=0; j<DCT_LEN; j++){
sum = 0;
for(k=0; k<DCT_LEN; k++) sum += dctTableA[i][k] * in[k][j];
tmp[i][j] = sum;
}
}
for(i=0; i<DCT_LEN; i++){
for(j=0; j<DCT_LEN; j++){
sum = 0;
for(k=0; k<DCT_LEN; k++) sum += dctTableA[j][k] * tmp[i][k];
out[i][j] = sum;
}
}
}
// 二维离散余弦逆变换
void dct_2Didct(float out[][DCT_LEN], float in[][DCT_LEN]){
float tmp[DCT_LEN][DCT_LEN];
float sum;
int i, j, k;
for(i=0; i<DCT_LEN; i++){
for(j=0; j<DCT_LEN; j++){
sum = 0;
for(k=0; k<DCT_LEN; k++) sum += dctTableA[k][i] * in[k][j];
tmp[i][j] = sum;
}
}
for(i=0; i<DCT_LEN; i++){
for(j=0; j<DCT_LEN; j++){
sum = 0;
for(k=0; k<DCT_LEN; k++) sum += dctTableA[k][j] * tmp[i][k];
out[i][j] = sum;
}
}
}
- 测试代码
//==============================================================================
//
// 作者: 王小康
// 描述: 一维二维DCT变换测试程序
// 日期: 2019-05-09
//
//==============================================================================
#include <stdio.h>
#include <stdlib.h>
#include "dct.h"
static void show1D(float *x){
int i;
for(i=0; i<DCT_LEN; i++){
printf("%15f", x[i]);
}
putchar('\n');
}
static void show2D(float x[][DCT_LEN]){
int i, j;
for(i=0; i<DCT_LEN; i++){
for(j=0; j<DCT_LEN; j++){
printf("%15f", x[i][j]);
}
putchar('\n');
}
putchar('\n');
}
int main(int argc, char *argv[]){
float data1[DCT_LEN][DCT_LEN] = {
{ 89, 101, 114, 125, 126, 115, 105, 96},
{ 97, 115, 131, 147, 149, 135, 123, 113},
{114, 134, 159, 178, 175, 164, 149, 137},
{121, 143, 177, 196, 201, 189, 165, 150},
{119, 141, 175, 201, 207, 186, 162, 144},
{107, 130, 165, 189, 192, 171, 144, 125},
{ 97, 119, 149, 171, 172, 145, 117, 96},
{ 88, 107, 136, 156, 155, 129, 97, 75}
// {200, 200, 200, 200, 200, 200, 200, 200},
// {240, 220, 190, 20, 20, 60, 30, 160},
// {200, 160, 180, 40, 40, 50, 50, 120},
// {220, 120, 190, 60, 60, 80, 70, 100},
// {190, 100, 100, 80, 80, 90, 99, 111},
// {100, 150, 120, 100, 70, 100, 120, 130},
// {111, 111, 121, 123, 110, 111, 134, 167},
// {120, 220, 222, 222, 210, 109, 190, 100}
};
float data2[DCT_LEN][DCT_LEN];
float data3[DCT_LEN][DCT_LEN];
dct_init();
printf("一维原始数据 ：");
show1D((float*)data1);
dct_1Dspectrum((float*)data2, (float*)data1);
printf("一维原始数据频谱：");
show1D((float*)data2);
dct_1Ddct((float*)data2, (float*)data1);
printf("一维DCT变换后 ：");
show1D((float*)data2);
dct_1Didct((float*)data3, (float*)data2);
printf("一维DCT逆变换后 ：");
show1D((float*)data3);
putchar('\n');
printf("二维原始数据：\n");
show2D(data1);
putchar('\n');
dct_2Ddct(data2, data1);
printf("二维DCT变换后：\n");
show2D(data2);
putchar('\n');
dct_2Didct(data3, data2);
printf("二维DCT逆变换后：\n");
show2D(data3);
putchar('\n');
return 0;
}
-
运行效果
