Computing Square Root of A Vector of Fixed-Point Numbers

August 9, 2011 Coded in ASM for the TI C64x
* =========================================================================== *
*                                                                             *
*  Compute square root of a Q.15 vector by 4th order polynomial fitting       *
*      y[i] = sqrt(x[i]), 0 <= x[i] < 1;                                      *
*                                                                             * 
*  C prototype:                                                               *
*      void DSP_vsqrt_q15(short* x, short* y, int N);                         *
*                                                                             * 
*  Performance:                                                               *
*      O(4*N) or 4 cycles per number (software pipelining enabled by -O2)     * 
*                                                                             * 
*  Error:                                                                     *
*      (-2^-15, 2^-15)                                                        *
*                                                                             *
* =========================================================================== *  

; chebfun, min max error        
; 0xFFF1024A = Q31(-0.0005)      
; 0x0046E0AB = Q31(0.0022)       
; 0xFE764EE1 = Q31(-0.0120)      
; 0x1277E288 = Q31(0.1443)       
; 0x6ED9EBA1 = Q31(0.8660)       

; polyfit, min sqr error   
; 0xFFF1203D = Q31(-0.0005)
; 0x004581E7 = Q31(0.0021) 
; 0xFE7645AF = Q31(-0.0120)
; 0x1278CF97 = Q31(0.1443) 
; 0x6ED9E687 = Q31(0.8660) 

SQRT_C4 .set 0xFFF1203D  
SQRT_C3 .set 0x004581E7  
SQRT_C2 .set 0xFE7645AF  
SQRT_C1 .set 0x1278CF97  
SQRT_C0 .set 0x6ED9E687  
SQRT_S  .set 0x5A82799A ;  Q31(0.7071)

        .sect ".text: _DSP_vsqrt_q15"
        .global _DSP_vsqrt_q15

_DSP_vsqrt_q15: .cproc  A_X, B_Y, A_n

        .no_mdep
        .rega A_C4, A_C3, A_C2, A_C1, A_C0, A_xx0, A_rnd, A_y0
        .rega A_y0c4, A_y0c3, A_y0c2, A_y0c1, A_x0c3, A_x0c2, A_x0c1, A_x0c0
        .rega A_S, A_e0, A_x0x, A_xm, A_y0l, A_y0h, A_y0s
        .regb B_C4, B_C3, B_C2, B_C1, B_C0, B_xx1, B_rnd, B_y1, B_y10
        .regb B_y1c4, B_y1c3, B_y1c2, B_y1c1, B_x1c3, B_x1c2, B_x1c1, B_x1c0
        .regb B_S, B_e1, B_x1x, B_xm, B_y1l, B_y1h, B_y1s, B_X
        .reg  B_i, C0, C1
        
            ADD         2,              A_X,            B_X
            MVK         0x1,            A_rnd
            SHL         A_rnd,          15,             A_rnd
            MV          A_rnd,          B_rnd
            
            MVKL        SQRT_C4,        A_C4
            MVKH        SQRT_C4,        A_C4
            MVKL        SQRT_C3,        B_C3
            MVKH        SQRT_C3,        B_C3

            MV          A_C4,           B_C4
            MV          B_C3,           A_C3
            MVKL        SQRT_C2,        A_C2
            MVKH        SQRT_C2,        A_C2
            MVKL        SQRT_C1,        B_C1
            MVKH        SQRT_C1,        B_C1

            MV          A_C2,           B_C2
            MV          B_C1,           A_C1
            MVKL        SQRT_C0,        A_C0
            MVKH        SQRT_C0,        A_C0
            MVKL        SQRT_S,         B_S
            MVKH        SQRT_S,         B_S
            
            MV          B_S,            A_S
            MV          A_C0,           B_C0

            MVKL        0x6000,         A_xm
            SHL         A_xm,           16,             A_xm
            MV          A_xm,           B_xm

            SHR         A_n,            1,              B_i
            SUB         B_i,            2,              B_i
                  
LOOP_vsqrt: .trip 8
            LDH.D1T1    *A_X++[2],      A_xx0
            LDH.D2T2    *B_X++[2],      B_xx1
            
            NORM.L1     A_xx0,          A_e0
            NORM.L2     B_xx1,          B_e1          
            SSHVL.M1    A_xx0,          A_e0,           A_x0x
            SSHVL.M2    B_xx1,          B_e1,           B_x1x

            SUB.D1      A_e0,           16,             A_e0
            SUB.D2      B_e1,           16,             B_e1
            AND.D1      0x1,            A_e0,           C0
            AND.D2      0x1,            B_e1,           C1
            SHR.S1      A_e0,           1,              A_e0
            SHR.S2      B_e1,           1,              B_e1

            SUB.D1      A_x0x,          A_xm,           A_x0x
            SUB.D2      B_x1x,          B_xm,           B_x1x
            SHL.S1      A_x0x,          2,              A_x0x
            SHL.S2      B_x1x,          2,              B_x1x

            MPYHIR.M1   A_x0x,          A_C4,           A_y0c4
            MPYHIR.M2   B_x1x,          B_C4,           B_y1c4
            SADD.L1     A_y0c4,         A_C3,           A_x0c3
            SADD.L2     B_y1c4,         B_C3,           B_x1c3
    
            MPYHIR.M1   A_x0x,          A_x0c3,         A_y0c3
            MPYHIR.M2   B_x1x,          B_x1c3,         B_y1c3
            SADD.L1     A_y0c3,         A_C2,           A_x0c2
            SADD.L2     B_y1c3,         B_C2,           B_x1c2

            MPYHIR.M1   A_x0x,          A_x0c2,         A_y0c2
            MPYHIR.M2   B_x1x,          B_x1c2,         B_y1c2
            SADD.L1     A_y0c2,         A_C1,           A_x0c1
            SADD.L2     B_y1c2,         B_C1,           B_x1c1

            MPYHIR.M1   A_x0x,          A_x0c1,         A_y0c1
            MPYHIR.M2   B_x1x,          B_x1c1,         B_y1c1
            SADD.L1     A_y0c1,         A_C0,           A_x0c0
            SADD.L2     B_y1c1,         B_C0,           B_x1c0

            ; A_S = B_S = 0x5A82799A ~= 0x5A820000 + 0x00008000
     [C0]   MPYHIR.M1   A_S,            A_x0c0,         A_y0h
     [C1]   MPYHIR.M2   B_S,            B_x1c0,         B_y1h
     [C0]   SHR         A_x0c0,         16,             A_y0l
     [C1]   SHR         B_x1c0,         16,             B_y1l 
     [C0]   SADD.L1     A_y0h,          A_y0l,          A_x0c0
     [C1]   SADD.L2     B_y1h,          B_y1l,          B_x1c0

            SHR.S1      A_x0c0,         A_e0,           A_y0s
            SHR.S2      B_x1c0,         B_e1,           B_y1s

            SADD.L1     A_y0s,          A_rnd,          A_y0
            SADD.L2     B_y1s,          B_rnd,          B_y1
            PACKH2.S2X  B_y1,           A_y0,           B_y10                       

            STW.D2T2    B_y10,          *B_Y++

            BDEC        LOOP_vsqrt,     B_i
            
            .endproc

First order RC Filter in the Digital Domain

August 7, 20111 comment Coded in Matlab
function [b,a] = rc_filter(R, C, Fs, filter_type)
% Returns equivalent IIR coefficients for an analog RC filter
%
% Usage:     [B,A] = RC_FILTER(r, c, fs, type);
%
%             R is the resistance value (in ohms)
%             C is the capacitance value (in farrads)
%             FS is the digital sample rate (in Hz)
%             type is a character string defining filter type
%                 Choices are: 'high' or 'low'
%
% Highpass filters have the analog structure:
%
%
%                  | | 
%     Vi  o--------| |----------+---------o  Vo
%                  | |          |
%                     C         | 
%                              --- 
%                              | |  R
%                              | | 
%                              ---   
%                               |
%                               |
%         o---------------------+---------o
%                              GND
%
%
% Lowpass filters have the analog structure:
%
%
%                  |-----| 
%     Vi  o--------|     |------+---------o  Vo
%                  |-----|      |
%                     R         | 
%                             ----- C
%                             ----- 
%                               | 
%                               |
%         o---------------------+---------o
%                              GND
%
% This function uses a pre-calculated equation for both of these circuits
% that only requires the resistance and capacitance value to get a true
% digital filter equivalent to a basic analog filter.
%
% The math behind these equations is based off the basic bilinear transform
% technique that can be found in many DSP textbooks.  The reference paper
% for this function was "Conversion of Analog to Digital Transfer 
% Functions" by C. Sidney Burrus, page 6.
%
% This is also the equivalent of a 1st order butterworth with a cuttoff
% frequency of Fc = 1/(2*pi*R*C);
%
% Author: sparafucile17 07/01/02
%

% Verify that cutoff of the analog filter is below Nyquist
if(  (1/(2*pi*R*C))  >  (Fs/2)  )
    error('This analog filter cannot be realized with this sample rate');
end

% Default to allpass if invalid type is selected
b = [ 1 0 ];
a = [ 1 0 ];

% Constants
RC = R * C;
T  = 1 / Fs;

% Analog Cutoff Fc
w = 1 / (RC);

% Prewarped coefficient for Bilinear transform
A = 1 / (tan((w*T) / 2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following equations were derived from
%
%            s
% T(s) =  -------
%          s + 1
%
%
% using Bilinear transform of
%
%             1          ( 1 - z^-1 )
% s -->  -----------  *  ------------
%         tan(w*T/2)     ( 1 + z^-1 )
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if(strcmp(filter_type,'high'))
    
    b(1) = (A)     / (1 + A);
    b(2) = -b(1);
    a(2) = (1 - A) / (1 + A);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following equations were derived from
%
%            1
% T(s) =  -------
%          s + 1
%
%
% using Bilinear transform of
%
%             1          ( 1 - z^-1 )
% s -->  -----------  *  ------------
%         tan(w*T/2)     ( 1 + z^-1 )
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if(strcmp(filter_type,'low'))
    
    b(1) = (1)     / (1 + A);
    b(2) = b(1);
    a(2) = (1 - A) / (1 + A);

end

Exponential Audio Unmute

August 7, 2011 Coded in Matlab
function [y] = signal_unmute(x, index, duration, Fs)
% Return the input vector with a unmute that occurs at a specific
% index and with an exponential ramp-up to reduce "pop" sounds.
% The output will start off at -100dB gain and end at 0dB gain.
%
% Usage: y = SIGNAL_UNMUTE(x, index, duration, Fs);
%
%        X is your one-dimensional input array
%        INDEX is where in the input signal you want the unmute to begin
%        DURATION is how long (in seconds) to exponentially ramp up the
%           input signal. 100ms is recommended.
%        FS is the sample rate
%
% Example: 
%        You want to unmute your signal at 0.5sec with a duration of 100ms 
%        and with a 48k Fs. (signal is unmuted completely at 0.6sec)
%        y = signal_mute(x, 24000, 0.1, 48000);
%
% Author: sparafucile17 7/29/2003

% Input must have some length
if(length(x) == 1)
    error('ERROR: input signal must have more than one element');
end

% This function only supports one-dimensional arrays
if((size(x, 2) ~= 1) && (size(x, 1) ~= 1))
    error('ERROR: Input must be one-dimensional');
end

% Make sure there are enough samples to complete the mute
if(length(x) < (index + duration*Fs))
    error(['There are not enough samples in X to complete the unmute.  '  ...
           'Either change the mute duration or move the index back.' ]);
end

% Flip vector (temporarily)
if((size(x, 2) ~= 1))
    x = x';
    flip_back = true;
else
    flip_back = false;
end

% Calculate exponential coefficient
dB_atten  = -100; %What do we consider "muted" to be in decibels
decayrate = -dB_atten / duration;
coeff = 1.0 - 10^(-decayrate/(20.0*Fs));

% Build the Gain array
gain = [zeros(index, 1); ones((length(x)-index), 1);];
b = [ coeff  0];
a = [ 1 (-1*(1-coeff))];
gain  = filter(b,a,gain);

% Apply Mute (gain) to the input signal
y = gain .* x;

% Flip the vector (if required)
if(flip_back == true);
    y = y';
end

Exponential Audio Mute

August 7, 2011 Coded in Matlab
function [y] = signal_mute(x, index, duration, Fs)
% Return the input vector with a mute that occurs at a specific
% index and with an exponential ramp-down to reduce "pop" sounds.
% The output will start off at 0dB gain and end at -100dB gain.
%
% Usage: y = SIGNAL_MUTE(x, index, duration, Fs);
%
%        X is your one-dimensional input array
%        INDEX is where in the input signal you want the mute to begin
%        DURATION is how long (in seconds) to exponentially ramp down the
%           input signal. 100ms is recommended.
%        FS is the sample rate
%
% Example: 
%        You want to mute your signal at 0.5sec with a duration of 100ms 
%        and with a 48k Fs. (mute complete at 0.6sec)
%        y = signal_mute(x, 24000, 0.1, 48000);
%
% Author: sparafucile17 7/29/2003

% Input must have some length
if(length(x) == 1)
    error('ERROR: input signal must have more than one element');
end

% This function only supports one-dimensional arrays
if((size(x, 2) ~= 1) && (size(x, 1) ~= 1))
    error('ERROR: Input must be one-dimensional');
end

% Make sure there are enough samples to complete the mute
if(length(x) < (index + duration*Fs))
    error(['There are not enough samples in X to complete the mute.  '  ...
           'Either change the mute duration or move the index back.' ]);
end

% Flip vector (temporarily)
if((size(x, 2) ~= 1))
    x = x';
    flip_back = true;
else
    flip_back = false;
end

% Calculate exponential coefficient
dB_atten  = -100; %How much attenuation will be at time: index + duration
decayrate = -dB_atten / duration;
coeff = 1.0 - 10^(-decayrate/(20.0*Fs));

% Build the Gain array
gain = [ones(index, 1); zeros((length(x)-index), 1);];
b = [ coeff  0];
a = [ 1 (-1*(1-coeff))];
gain  = filter(b,a,gain, [1]);

% Remove minor overshoot at the beginning of gain vector
gain(find(gain > 1)) = 1;

% Apply Mute (gain) to the input signal
y = gain .* x;

% Flip the vector (if required)
if(flip_back == true);
    y = y';
end

Fast MDCT/IMDCT Based on Forward FFT

August 6, 2011 Coded in C
/******** begin of mdct.h ******** */
#ifndef __MDCT_H
#define __MDCT_H

#include <fftw3.h>

#ifdef __cplusplus
extern "C" {
#endif

#ifdef SINGLE_PRECISION
typedef float         FLOAT;
typedef fftwf_complex FFTW_COMPLEX;
typedef fftwf_plan    FFTW_PLAN;
#else // DOUBLE_PRECISION
typedef double        FLOAT;
typedef fftw_complex  FFTW_COMPLEX;
typedef fftw_plan     FFTW_PLAN;
#endif  // SINGLE_PRECISION

typedef struct {
    int           M;            // MDCT spectrum size (number of bins)
    FLOAT*        twiddle;      // twiddle factor
    FFTW_COMPLEX* fft_in;       // fft workspace, input
    FFTW_COMPLEX* fft_out;      // fft workspace, output
    FFTW_PLAN     fft_plan;     // fft configuration
} mdct_plan; 

mdct_plan* mdct_init(int M);    // MDCT spectrum size (number of bins)

void mdct_free(mdct_plan* m_plan);

void mdct(FLOAT* mdct_line, FLOAT* time_signal, mdct_plan* m_plan);

void imdct(FLOAT* time_signal, FLOAT* mdct_line, mdct_plan* m_plan);

#ifdef __cplusplus
}
#endif

#endif // __MDCT_H
/******** end of mdct.h ******** */

/******** begin of mdct.c ******** */
#ifdef SINGLE_PRECISION

#define FFTW_MALLOC   fftwf_malloc
#define FFTW_FREE     fftwf_free
#define FFTW_PLAN_1D  fftwf_plan_dft_1d
#define FFTW_DESTROY  fftwf_destroy_plan 
#define FFTW_EXECUTE  fftwf_execute

#else // DOUBLE_PRECISION

#define FFTW_MALLOC   fftw_malloc
#define FFTW_FREE     fftw_free
#define FFTW_PLAN_1D  fftw_plan_dft_1d
#define FFTW_DESTROY  fftw_destroy_plan 
#define FFTW_EXECUTE  fftw_execute

#endif // SINGLE_PRECISION

void mdct_free(mdct_plan* m_plan)
{
    if(m_plan) 
    {
        FFTW_DESTROY(m_plan->fft_plan);
        FFTW_FREE(m_plan->fft_in);
        FFTW_FREE(m_plan->fft_out);

        if(m_plan->twiddle)
            free(m_plan->twiddle);

        free(m_plan);
    }
}

#define MDCT_CLEAUP(msg, ...) \
    {fprintf(stderr, msg", %s(), %s:%d \n", \
            __VA_ARGS__, __func__, __FILE__, __LINE__); \
     mdct_free(m_plan); return NULL;}

mdct_plan* mdct_init(int M)
{
    int        n;
    FLOAT      alpha, omega, scale;
    mdct_plan* m_plan = NULL;

    if(0x00 != (M & 0x01)) 
        MDCT_CLEAUP(" Expect an even number of MDCT coeffs, but meet %d", M);

    m_plan = (mdct_plan*) malloc(sizeof(mdct_plan));
    if(NULL == m_plan)
        MDCT_CLEAUP(" malloc error: %s", "m_plan");
    memset(m_plan, 0, sizeof(m_plan[0]));

    m_plan->M = M;

    m_plan->twiddle = (FLOAT*) malloc(sizeof(FLOAT) * M);
    if(NULL == m_plan->twiddle)
        MDCT_CLEAUP(" malloc error: %s", "m_plan->twiddle");
    alpha = M_PI / (8.f * M);
    omega = M_PI / M;
    scale = sqrt(sqrt(2.f / M));    
    for(n = 0; n < (M >> 1); n++)
    {    
        m_plan->twiddle[2*n+0] = (FLOAT) (scale * cos(omega * n + alpha));
        m_plan->twiddle[2*n+1] = (FLOAT) (scale * sin(omega * n + alpha));
    }    
	
    m_plan->fft_in   
        = (FFTW_COMPLEX*) FFTW_MALLOC(sizeof(FFTW_COMPLEX) * M >> 1);    
    if(NULL == m_plan->fft_in)
        MDCT_CLEAUP(" malloc error: %s", "m_plan->fft_in");

    m_plan->fft_out  
        = (FFTW_COMPLEX*) FFTW_MALLOC(sizeof(FFTW_COMPLEX) * M >> 1);    
    if(NULL == m_plan->fft_out)
        MDCT_CLEAUP(" malloc error: %s", "m_plan->fft_out");

    m_plan->fft_plan = FFTW_PLAN_1D(M >> 1, 
                                    m_plan->fft_in, 
                                    m_plan->fft_out,    
                                    FFTW_FORWARD,
                                    FFTW_MEASURE);
    if(NULL == m_plan->fft_plan)
        MDCT_CLEAUP(" malloc error: %s", "m_plan->fft_plan");

    return m_plan;
}

void mdct(FLOAT* mdct_line, FLOAT* time_signal, mdct_plan* m_plan)
{
    FLOAT *xr, *xi, r0, i0;
    FLOAT *cos_tw, *sin_tw, c, s;
    int    M, M2, M32, M52, n;

    M   = m_plan->M;
    M2  = M >> 1;
    M32 = 3 * M2;
    M52 = 5 * M2;

    cos_tw = m_plan->twiddle;
    sin_tw = cos_tw + 1; 
    
    /* odd/even folding and pre-twiddle */
    xr = (FLOAT*) m_plan->fft_in;
    xi = xr + 1;
    for(n = 0; n < M2; n += 2) 
    {
        r0 = time_signal[M32-1-n] + time_signal[M32+n];    
        i0 = time_signal[M2+n]    - time_signal[M2-1-n];    
        
        c = cos_tw[n];
        s = sin_tw[n];

        xr[n] = r0 * c + i0 * s;
        xi[n] = i0 * c - r0 * s;
    }

    for(; n < M; n += 2) 
    {
        r0 = time_signal[M32-1-n] - time_signal[-M2+n];    
        i0 = time_signal[M2+n]    + time_signal[M52-1-n];    
        
        c = cos_tw[n];
        s = sin_tw[n];

        xr[n] = r0 * c + i0 * s;
        xi[n] = i0 * c - r0 * s;
    }

    /* complex FFT of size M/2 */
    FFTW_EXECUTE(m_plan->fft_plan);

    /* post-twiddle */
    xr = (FLOAT*) m_plan->fft_out;
    xi = xr + 1;
    for(n = 0; n < M; n += 2)
    {
        r0 = xr[n];
        i0 = xi[n];
        
        c = cos_tw[n];
        s = sin_tw[n];    

        mdct_line[n]     = - r0 * c - i0 * s;
        mdct_line[M-1-n] = - r0 * s + i0 * c;
    }
}

void imdct(FLOAT* time_signal, FLOAT* mdct_line, mdct_plan* m_plan)
{
    FLOAT *xr, *xi, r0, i0, r1, i1;
    FLOAT *cos_tw, *sin_tw, c, s;
    int    M, M2, M32, M52, n;

    M   = m_plan->M;
    M2  = M >> 1;
    M32 = 3 * M2;
    M52 = 5 * M2;

    cos_tw = m_plan->twiddle;
    sin_tw = cos_tw + 1; 
    
    /* pre-twiddle */
    xr = (FLOAT*) m_plan->fft_in;
    xi = xr + 1;
    for(n = 0; n < M; n += 2)
    {
        r0 =  mdct_line[n];	
        i0 =  mdct_line[M-1-n];
        
        c = cos_tw[n];
        s = sin_tw[n];    
        
        xr[n] = -i0 * s - r0 * c;
        xi[n] = -i0 * c + r0 * s;
    } 
    
    /* complex FFT of size M/2 */
    FFTW_EXECUTE(m_plan->fft_plan);

    /* odd/even expanding and post-twiddle */
    xr = (FLOAT*) m_plan->fft_out;
    xi = xr + 1;
    for(n = 0; n < M2; n += 2) 
    {
        r0 = xr[n]; 
        i0 = xi[n];    
        
        c = cos_tw[n]; 
        s = sin_tw[n];

        r1 = r0 * c + i0 * s;
        i1 = r0 * s - i0 * c;

        time_signal[M32-1-n] =  r1;
        time_signal[M32+n]   =  r1;
        time_signal[M2+n]    =  i1;
        time_signal[M2-1-n]  = -i1;
    }

    for(; n < M; n += 2) 
    {
        r0 = xr[n]; 
        i0 = xi[n];    
        
        c = cos_tw[n]; 
        s = sin_tw[n];
        
        r1 = r0 * c + i0 * s;
        i1 = r0 * s - i0 * c;
        
        time_signal[M32-1-n] =  r1;
        time_signal[-M2+n]   = -r1;
        time_signal[M2+n]    =  i1;
        time_signal[M52-1-n] =  i1;
    }
}
/******** end of mdct.c ******** */

/******** begin of mdct_test.c ******** */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>
#include <time.h>
#include "mdct.h"

int main(int argc, char* argv[])
{
    int        M, r, i;
    FLOAT*     time = NULL;
    FLOAT*     freq = NULL;
    mdct_plan* m_plan = NULL;
    char*      precision = NULL;
    struct timeval t0, t1;
    long long elps;

    if(3 != argc)
    {
        fprintf(stderr, " Usage: %s <MDCT_SPECTRUM_SIZE> <run_times> \n", argv[0]);
        return -1;
    }

    sscanf(argv[1], "%d", &M);
    sscanf(argv[2], "%d", &r);    
    if(NULL == (m_plan = mdct_init(M)))
        return -1;
    if(NULL == (time = (FLOAT*) malloc(2 * M * sizeof(FLOAT))))
        return -1;
    if(NULL == (freq = (FLOAT*) malloc(M * sizeof(FLOAT))))
        return -1;

    for(i = 0; i < 2 * M; i++)
        time[i] = 2.f * rand() / RAND_MAX - 1.f;        
    for(i = 0; i < M; i++)
        freq[i] = 2.f * rand() / RAND_MAX - 1.f;        

    precision = (sizeof(float) == sizeof(FLOAT))? 
                "single precision" : "double precision";

#if 1
    gettimeofday(&t0, NULL);
    for(i = 0; i < r; i++)
        mdct(freq, time, m_plan);
    gettimeofday(&t1, NULL);

    elps = (t1.tv_sec - t0.tv_sec) * 1000000 + (t1.tv_usec - t0.tv_usec);    
    fprintf(stdout, "MDCT size of %d, %s, running %d times, average %.3f ms\n", 
            M, precision, r, (FLOAT) elps / r / 1000.f); 
#endif // 0

#if 1    
    gettimeofday(&t0, NULL);
    for(i = 0; i < r; i++)
        imdct(time, freq, m_plan);
    gettimeofday(&t1, NULL);

    elps = (t1.tv_sec - t0.tv_sec) * 1000000 + (t1.tv_usec - t0.tv_usec);    
    fprintf(stdout, "IMDCT size of %d, %s, running %d times, average %.3f ms\n", 
            M, precision, r, (FLOAT) elps / r / 1000.f); 
#endif //0

#if 0
    for(i = 0; i < 2 * M; i++)
        fprintf(stdout, "%f    ", time[i]);
    fprintf(stdout, "\n");

    for(i = 0; i < M; i++)
        fprintf(stdout, "%f    ", freq[i]);
    fprintf(stdout, "\n");
#endif // 0

    free(time);
    free(freq);
    mdct_free(m_plan);    

    return 0;
}
/******** end of mdct_test.c ******** */

Image denoising: Threshold calculation for vishushrink method

Senthilkumar July 30, 2011 Coded in Matlab
%function to calculate the threshold using Visushrink denoising method
function T = Visu_threshold(X)
    [m,n] = size(X);
    M = m*n;
    T = sqrt(2*log(M));

Image denoising: using vishushrink method

Senthilkumar July 30, 2011 Coded in Matlab
function [soft_X1,SOFT_PSNR] = Vishu_soft(X,Y)
%function used to denoise a noisy image using vishuShrink method
%One -level decomposition
[CA,CH,CV,CD] = dwt2(Y,'haar');
%Call the function to calculate the threshold
T1 = Visu_threshold(CD);
% Call the function to perfom soft shrinkage
de_CH = soft(CH,T1);
de_CV = soft(CV,T1);
de_CD = soft(CD,T1);
%
%Two -level decomposition
[CA1,CH1,CV1,CD1] = dwt2(CA,'haar');
%Call the function to calculate the threshold
T2 = Visu_threshold(CD1);
% Call the function to perfom soft shrinkage
de_CH1 = soft(CH1,T2);
de_CV1 = soft(CV1,T2);
de_CD1 = soft(CD1,T2);
% % CA1 = soft1(CA1,T2);
%
%
%Reconstruction for soft shrinkage
X2 = idwt2(CA1,de_CH1,de_CV1,de_CD1,'haar');
X1 = idwt2(X2,de_CH,de_CV,de_CD,'haar');

SOFT_PSNR = PSNR(X,X1);
soft_X1 = uint8(X1);

Image denoise: Threshold for sureshrink method

Senthilkumar July 30, 20111 comment Coded in Matlab
function thr = sureshrink(CD,T)
%function used to calculate threshold using sureshrink method
CD = CD(:)';
n = length(CD);
sx2 = sort(abs(CD)-T).^2;  % sorting in descending order
b  = cumsum(sx2);           %cumulative sum
risks = (n-(2*(1:n))+b)/n;
[risk,best] = min(risks);
thr = sqrt(sx2(best));

Image denoise: denoising using soft threshold

Senthilkumar July 30, 2011 Coded in Matlab
function y = soft(x,T)
%function used to denoise a noisy image with given soft threshold
%x = noisy image
%T = threshold
y = max(abs(x) - T, 0);
y = y./(y+T) .* x;

Image denoising: Peak signal to Noise ratio calculation

Senthilkumar July 30, 20111 comment Coded in Matlab
function A = PSNR(G,H)
%G = original image
%H =  denoised image
  error = G - H;
  decibels = 20*log10((255*255)/(sqrt(mean(mean(error.^2)))));
  disp(sprintf('PSNR = %5.2f db',decibels)) 
  A = decibels;