请高手帮我解决这个离散余弦变换(DCT)的问题。高分相送。

void dct(FP  *s,
         int N,
int opt)
{
int n, k;
FP alpha;
FP pin, sqrtn, sqrt2n;
FP *x, *X;
FP *sp, *xp, *Xp;

pin    = 0.5*PI/N;
sqrtn  = sqrt(1.0/N);
sqrt2n = sqrt(2.0/N);

if (opt==DCT) { // Forward DCT
/**** SOURCE CODE MISSING !!!! ****/
// Calculate the real-valued DCT       
 

}
else { // inverse DCT
/**** SOURCE CODE MISSING !!!! ****/
// Calculate the real-valued IDCT

}

return;
}

其中*s是既是输入，也是完成离散余弦变换计算后的输出，N是要计算点的个数。alpha是系数。当k等于0时候，alpha=sqrtn,k不等于0的时候k=sqrt2n.
下面有DCT和IDCT的方程
http://documents.wolfram.com/applications/digitalimage/UsersGuide/8.4.html
请参考。其中的n,k和程序中的都一样。你不用考虑opt,当opt=1的时候是DCT。补全程序就可以了。

5 个解决方案

#1

这个网站是国外的，俺这边进不去，参考一下下面的吧

/************************************DCT**************************************/
//  Set up the cosine coefficients c0..c7.           
#define C0 0xB505
#define C1 0x2C62
#define C2 0x29CF
#define C3 0x25A0
#define C5 0x1924
#define C6 0x1151
#define C7 0x08D4

void fdct(short *dct_data)//, unsigned num_fdcts)
{                                                              
    int   g0, g1, h0, h1,                                        
          p0, p1;               // Even-half intermediate. /
    short r0, r1;               // Even-half intermediate. /
    int   P0, P1, R0, R1;       // Even-half intermediate. /
    short g2, g3, h2, h3;       // Odd-half intermediate.  /
    short q0a,s0a,q0, q1,                                       
          s0, s1;               // Odd-half intermediate.  /
    short Q0, Q1, S0, S1;       // Odd-half intermediate.  /
    int   F0, F1, F2, F3,                                        
          F4, F5, F6, F7;       // Freq. domain results.   /
    unsigned i, j;               
    short *dataptr;                                        
                                                                 
    //  Perform Vert 1-D FDCT on columns within each block. 
dataptr = dct_data;                                       
    for (j = 0; j < 8; j++)                                  
    {                                                        
        //  Stage 1:  Separate into even and odd halves./
        g0 = dataptr[ 0] + dataptr[56]; h2 = dataptr[ 0] - dataptr[56];            
        g1 = dataptr[ 8] + dataptr[48]; h3 = dataptr[ 8] - dataptr[48];            
        h1 = dataptr[16] + dataptr[40]; g3 = dataptr[16] - dataptr[40];            
        h0 = dataptr[24] + dataptr[32]; g2 = dataptr[24] - dataptr[32];            
                                                             
        //  Stage 2
        p0 = g0 + h0;               r0 = g0 - h0;            
        p1 = g1 + h1;               r1 = g1 - h1;            
        q1 = g2;                    s1 = h2;                 
                                                             
        s0a= h3 + g3;               q0a= h3 - g3;            
        s0 = (s0a * C0 + 0x7FFF) >> 16;                      
        q0 = (q0a * C0 + 0x7FFF) >> 16;                      
                                                             
        //  Stage 3 
        P0 = p0 + p1;               P1 = p0 - p1;            
        R1 = C6 * r1 + C2 * r0;     R0 = C6 * r0 - C2 * r1;  
                                                             
        Q1 = q1 + q0;               Q0 = q1 - q0;            
        S1 = s1 + s0;               S0 = s1 - s0;            
                                                             
        //  Stage 4
        F0 = P0;                    F4 = P1;                 
        F2 = R1;                    F6 = R0;                 
                                                             
        F1 = C7 * Q1 + C1 * S1;     F7 = C7 * S1 - C1 * Q1;  
        F5 = C3 * Q0 + C5 * S0;     F3 = C3 * S0 - C5 * Q0;  
                                                             
        //  Store the frequency domain results. 
        dataptr[ 0] = F0;                                 
        dataptr[ 8] = F1 >> 13;                           
        dataptr[16] = F2 >> 13;                           
        dataptr[24] = F3 >> 13;                           
        dataptr[32] = F4;                                 
        dataptr[40] = F5 >> 13;                           
        dataptr[48] = F6 >> 13;                           
        dataptr[56] = F7 >> 13;                           
                                                             
        dataptr++;                                        
    }                                                        
                                                                
    //  Perform Horizontal 1-D FDCT on each 8x8 block.    
    dataptr = dct_data;                                       
    for (i = 0; i < 8 ; i++)                          
    {                                                            
        //  Load the spatial-domain samples.               
        g0 = dataptr[0] + dataptr[7]; h2 = dataptr[0] - dataptr[7];            
        g1 = dataptr[1] + dataptr[6]; h3 = dataptr[1] - dataptr[6];            
        h1 = dataptr[2] + dataptr[5]; g3 = dataptr[2] - dataptr[5];            
        h0 = dataptr[3] + dataptr[4]; g2 = dataptr[3] - dataptr[4];            
                                                                 
        //  Stage 2                                       
        p0 = g0 + h0;               r0 = g0 - h0;                
        p1 = g1 + h1;               r1 = g1 - h1;                
        q1 = g2;                    s1 = h2;                     
                                                                 
        s0a= h3 + g3;               q0a= h3 - g3;                
        q0 = (q0a * C0 + 0x7FFF) >> 16;                          
        s0 = (s0a * C0 + 0x7FFF) >> 16;                          
                                                                 
        //  Stage 3                                        
        P0 = p0 + p1;               P1 = p0 - p1;                
        R1 = C6 * r1 + C2 * r0;     R0 = C6 * r0 - C2 * r1;      
                                                                 
        Q1 = q1 + q0;               Q0 = q1 - q0;                
        S1 = s1 + s0;               S0 = s1 - s0;                
                                                                 
        //  Stage 4                                        
        F0 = P0;                    F4 = P1;                     
        F2 = R1;                    F6 = R0;                     
                                                                 
        F1 = C7 * Q1 + C1 * S1;     F7 = C7 * S1 - C1 * Q1;      
        F5 = C3 * Q0 + C5 * S0;     F3 = C3 * S0 - C5 * Q0;      
                                                                 
        //  Round and truncate values.                     
        //  Note: F0 and F4 have different rounding since no
        //  MPYs have been applied to either term.  Also, F0's  
        //  rounding is slightly different to offset the  
        //  truncation effects from the horizontal pass (which 
        //  does not round).                              

dataptr[0] = (F0 + 0x0006) >>  3;                                      
        dataptr[1] = (F1 + 0x7FFF) >> 16;                                 
        dataptr[2] = (F2 + 0x7FFF) >> 16;                                   
        dataptr[3] = (F3 + 0x7FFF) >> 16;                                    
        dataptr[4] = (F4 + 0x0004) >>  3;                                 
        dataptr[5] = (F5 + 0x7FFF) >> 16;                                  
        dataptr[6] = (F6 + 0x7FFF) >> 16;                                   
        dataptr[7] = (F7 + 0x7FFF) >> 16;  
        //  Update pointer to next FDCT row. 
        dataptr += 8;                                         
    }                                                            
}           

#2

///////////////////////////////////////////////////////////////////////
void Initialize_Fast_IDCT()
{
short i;

iclp = iclip+512;
for (i= -512; i<512; i++)
iclp[i] = (i<-256) ? -256 : ((i>255) ? 255 : i);
}
////////////////////////////////////////////////////////////////////////
void idctrow(int * blk)
{
int x0, x1, x2, x3, x4, x5, x6, x7, x8;
//intcut
if (!((x1 = blk[4]<<11) | (x2 = blk[6]) | (x3 = blk[2]) |
(x4 = blk[1]) | (x5 = blk[7]) | (x6 = blk[5]) | (x7 = blk[3])))
{
blk[0]=blk[1]=blk[2]=blk[3]=blk[4]=blk[5]=blk[6]=blk[7]=blk[0]<<3;
return;
}
x0 = (blk[0]<<11) + 128; // for proper rounding in the fourth stage 
//first stage
x8 = W7*(x4+x5);
x4 = x8 + (W1-W7)*x4;
x5 = x8 - (W1+W7)*x5;
x8 = W3*(x6+x7);
x6 = x8 - (W3-W5)*x6;
x7 = x8 - (W3+W5)*x7;
//second stage
x8 = x0 + x1;
x0 -= x1;
x1 = W6*(x3+x2);
x2 = x1 - (W2+W6)*x2;
x3 = x1 + (W2-W6)*x3;
x1 = x4 + x6;
x4 -= x6;
x6 = x5 + x7;
x5 -= x7;
//third stage
x7 = x8 + x3;
x8 -= x3;
x3 = x0 + x2;
x0 -= x2;
x2 = (181*(x4+x5)+128)>>8;
x4 = (181*(x4-x5)+128)>>8;
//fourth stage
blk[0] = (x7+x1)>>8;
blk[1] = (x3+x2)>>8;
blk[2] = (x0+x4)>>8;
blk[3] = (x8+x6)>>8;
blk[4] = (x8-x6)>>8;
blk[5] = (x0-x4)>>8;
blk[6] = (x3-x2)>>8;
blk[7] = (x7-x1)>>8;
}
//////////////////////////////////////////////////////////////////////////////
void idctcol(int * blk)
{
int x0, x1, x2, x3, x4, x5, x6, x7, x8;
//intcut
if (!((x1 = (blk[8*4]<<8)) | (x2 = blk[8*6]) | (x3 = blk[8*2]) |
(x4 = blk[8*1]) | (x5 = blk[8*7]) | (x6 = blk[8*5]) | (x7 = blk[8*3])))
{
blk[8*0]=blk[8*1]=blk[8*2]=blk[8*3]=blk[8*4]=blk[8*5]
=blk[8*6]=blk[8*7]=iclp[(blk[8*0]+32)>>6];
return;
}
x0 = (blk[8*0]<<8) + 8192;
//first stage
x8 = W7*(x4+x5) + 4;
x4 = (x8+(W1-W7)*x4)>>3;
x5 = (x8-(W1+W7)*x5)>>3;
x8 = W3*(x6+x7) + 4;
x6 = (x8-(W3-W5)*x6)>>3;
x7 = (x8-(W3+W5)*x7)>>3;
//second stage
x8 = x0 + x1;
x0 -= x1;
x1 = W6*(x3+x2) + 4;
x2 = (x1-(W2+W6)*x2)>>3;
x3 = (x1+(W2-W6)*x3)>>3;
x1 = x4 + x6;
x4 -= x6;
x6 = x5 + x7;
x5 -= x7;
//third stage
x7 = x8 + x3;
x8 -= x3;
x3 = x0 + x2;
x0 -= x2;
x2 = (181*(x4+x5)+128)>>8;
x4 = (181*(x4-x5)+128)>>8;
//fourth stage
blk[8*0] = iclp[(x7+x1)>>14];
blk[8*1] = iclp[(x3+x2)>>14];
blk[8*2] = iclp[(x0+x4)>>14];
blk[8*3] = iclp[(x8+x6)>>14];
blk[8*4] = iclp[(x8-x6)>>14];
blk[8*5] = iclp[(x0-x4)>>14];
blk[8*6] = iclp[(x3-x2)>>14];
blk[8*7] = iclp[(x7-x1)>>14];
}

void Fast_IDCT(int * block)
{
short i;

for (i=0; i<8; i++)
idctrow(block+8*i);

for (i=0; i<8; i++)
idctcol(block+i);
}

#3

是作业吧
帮你顶

#4

天，我已经忘记数学了

#5

看不懂，顶～