DCT变换(离散余弦变换)

• 离散余弦变换原理
  一句话概括核心原理：频率不同的正弦波相乘，其黎曼积分总是为零。下图是网上找来的离散余弦变换公式：
  
• 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;
}

• 运行效果