mirdwt_r.c
上传用户:speoil
上传日期:2022-06-23
资源大小:224k
文件大小:7k
源码类别:

波变换

开发平台:

Matlab

  1. /*
  2. File Name: MIRDWT.c
  3. Last Modification Date: 06/14/95 16:22:45
  4. Current Version: MIRDWT.c 2.4
  5. File Creation Date: Wed Oct 12 08:44:43 1994
  6. Author: Markus Lang  <lang@jazz.rice.edu>
  7. Copyright (c) 2000 RICE UNIVERSITY. All rights reserved.
  8. Created by Markus Lang, Department of ECE, Rice University. 
  9. This software is distributed and licensed to you on a non-exclusive 
  10. basis, free-of-charge. Redistribution and use in source and binary forms, 
  11. with or without modification, are permitted provided that the following 
  12. conditions are met:
  13. 1. Redistribution of source code must retain the above copyright notice, 
  14.    this list of conditions and the following disclaimer.
  15. 2. Redistribution in binary form must reproduce the above copyright notice, 
  16.    this list of conditions and the following disclaimer in the 
  17.    documentation and/or other materials provided with the distribution.
  18. 3. All advertising materials mentioning features or use of this software 
  19.    must display the following acknowledgment: This product includes 
  20.    software developed by Rice University, Houston, Texas and its contributors.
  21. 4. Neither the name of the University nor the names of its contributors 
  22.    may be used to endorse or promote products derived from this software 
  23.    without specific prior written permission.
  24. THIS SOFTWARE IS PROVIDED BY WILLIAM MARSH RICE UNIVERSITY, HOUSTON, TEXAS, 
  25. AND CONTRIBUTORS AS IS AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, 
  26. BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS 
  27. FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL RICE UNIVERSITY 
  28. OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 
  29. EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 
  30. PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; 
  31. OR BUSINESS INTERRUPTIONS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, 
  32. WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR 
  33. OTHERWISE), PRODUCT LIABILITY, OR OTHERWISE ARISING IN ANY WAY OUT OF THE 
  34. USE OF THIS SOFTWARE,  EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  35. For information on commercial licenses, contact Rice University's Office of 
  36. Technology Transfer at techtran@rice.edu or (713) 348-6173
  37. Change History: Fixed the code such that 1D vectors passed to it can be in
  38.                 either passed as a row or column vector. Also took care of 
  39. the code such that it will compile with both under standard
  40. C compilers as well as for ANSI C compilers
  41. Jan Erik Odegard <odegard@ece.rice.edu> Wed Jun 14 1995
  42.                 Fix minor bug to allow maximum number of levels
  43. MATLAB description:
  44. %function x = mirdwt(yl,yh,h,L);
  45. % function computes the inverse redundant discrete wavelet transform y for a
  46. % 1D or  2D input signal. redundant means here that the subsampling after
  47. % each stage of the forward transform has been omitted. yl contains the
  48. % lowpass and yl the highpass components as computed, e.g., by mrdwt. In
  49. % case of a 2D signal the ordering in yh is [lh hl hh lh hl ... ] (first
  50. % letter refers to row, second to column filtering).  
  51. %
  52. %    Input:
  53. %       yl   : lowpass component
  54. %       yh   : highpass components
  55. %       h    : scaling filter
  56. %       L    : number of levels. in case of a 1D signal length(yl) must be
  57. %              divisible by 2^L; in case of a 2D signal the row and the
  58. %              column dimension must be divisible by 2^L.
  59. %   
  60. %    Output:
  61. % x    : finite length 1D or 2D signal
  62. %
  63. % see also: mdwt, midwt, mrdwt
  64. */
  65. #include <math.h>
  66. #include <stdio.h>
  67. #define max(a, b) ((a) > (b) ? (a) : (b))
  68. #define mat(a, i, j) (*(a + (m*(j)+i)))  /* macro for matrix indices */
  69. #ifdef __STDC__
  70. MIRDWT(double *x, int m, int n, double *h, int lh, int L,
  71.        double *yl, double *yh)
  72. #else
  73. MIRDWT(x, m, n, h, lh, L, yl, yh)
  74. double *x, *h, *yl, *yh;
  75. int m, n, lh, L;
  76. #endif
  77. {
  78.   double  *g0, *g1, *ydummyll, *ydummylh, *ydummyhl;
  79.   double *ydummyhh, *xdummyl , *xdummyh, *xh;
  80.   long i, j;
  81.   int actual_L, actual_m, actual_n, c_o_a, ir, n_c, n_cb, n_c_o, lhm1;
  82.   int ic, n_r, n_rb, n_r_o, c_o_a_p2n, sample_f;
  83.   xh = (double *)mxCalloc(m*n,sizeof(double));
  84.   xdummyl = (double *)mxCalloc(max(m,n),sizeof(double));
  85.   xdummyh = (double *)mxCalloc(max(m,n),sizeof(double));
  86.   ydummyll = (double *)mxCalloc(max(m,n)+lh-1,sizeof(double));
  87.   ydummylh = (double *)mxCalloc(max(m,n)+lh-1,sizeof(double));
  88.   ydummyhl = (double *)mxCalloc(max(m,n)+lh-1,sizeof(double));
  89.   ydummyhh = (double *)mxCalloc(max(m,n)+lh-1,sizeof(double));
  90.   g0 = (double *)mxCalloc(lh,sizeof(double));
  91.   g1 = (double *)mxCalloc(lh,sizeof(double));
  92.   
  93.   if (n==1){
  94.     n = m;
  95.     m = 1;
  96.   }
  97.   /* analysis lowpass and highpass */
  98.   for (i=0; i<lh; i++){
  99.     g0[i] = h[i]/2;
  100.     g1[i] = h[lh-i-1]/2;
  101.   }
  102.   for (i=1; i<=lh; i+=2)
  103.     g1[i] = -g1[i];
  104.   
  105.   lhm1 = lh - 1;
  106.   /* 2^L */
  107.   sample_f = 1;
  108.   for (i=1; i<L; i++)
  109.     sample_f = sample_f*2;
  110.   actual_m = m/sample_f;
  111.   actual_n = n/sample_f;
  112.   /* restore yl in x */
  113.   for (i=0;i<m*n;i++)
  114.     x[i] = yl[i];
  115.   
  116.   /* main loop */
  117.   for (actual_L=L; actual_L >= 1; actual_L--){
  118.     /* actual (level dependent) column offset */
  119.     if (m==1)
  120.       c_o_a = n*(actual_L-1);
  121.     else
  122.       c_o_a = 3*n*(actual_L-1);
  123.     c_o_a_p2n = c_o_a + 2*n;
  124.     
  125.     /* go by columns in case of a 2D signal*/
  126.     if (m>1){
  127.       n_rb = m/actual_m;                 /* # of row blocks per column */
  128.       for (ic=0; ic<n; ic++){            /* loop over column */
  129. for (n_r=0; n_r<n_rb; n_r++){    /* loop within one column */
  130.   /* store in dummy variables */
  131.   ir = -sample_f + n_r;
  132.   for (i=0; i<actual_m; i++){    
  133.     ir = ir + sample_f;
  134.     ydummyll[i+lhm1] = mat(x, ir, ic);  
  135.     ydummylh[i+lhm1] = mat(yh, ir, c_o_a+ic);  
  136.     ydummyhl[i+lhm1] = mat(yh, ir,c_o_a+n+ic);  
  137.     ydummyhh[i+lhm1] = mat(yh, ir, c_o_a_p2n+ic);   
  138.   }
  139.   /* perform filtering and adding: first LL/LH, then HL/HH */
  140.   bpconv(xdummyl, actual_m, g0, g1, lh, ydummyll, ydummylh); 
  141.   bpconv(xdummyh, actual_m, g0, g1, lh, ydummyhl, ydummyhh); 
  142.   /* store dummy variables in matrices */
  143.   ir = -sample_f + n_r;
  144.   for (i=0; i<actual_m; i++){    
  145.     ir = ir + sample_f;
  146.     mat(x, ir, ic) = xdummyl[i];  
  147.     mat(xh, ir, ic) = xdummyh[i];  
  148.   }
  149. }
  150.       }
  151.     }
  152.     
  153.     /* go by rows */
  154.     n_cb = n/actual_n;                 /* # of column blocks per row */
  155.     for (ir=0; ir<m; ir++){            /* loop over rows */
  156.       for (n_c=0; n_c<n_cb; n_c++){    /* loop within one row */      
  157. /* store in dummy variable */
  158. ic = -sample_f + n_c;
  159. for  (i=0; i<actual_n; i++){    
  160.   ic = ic + sample_f;
  161.   ydummyll[i+lhm1] = mat(x, ir, ic);  
  162.   if (m>1)
  163.     ydummyhh[i+lhm1] = mat(xh, ir, ic);  
  164.   else
  165.     ydummyhh[i+lhm1] = mat(yh, ir, c_o_a+ic);  
  166. /* perform filtering lowpass/highpass */
  167. bpconv(xdummyl, actual_n, g0, g1, lh, ydummyll, ydummyhh); 
  168. /* restore dummy variables in matrices */
  169. ic = -sample_f + n_c;
  170. for (i=0; i<actual_n; i++){    
  171.   ic = ic + sample_f;
  172.   mat(x, ir, ic) = xdummyl[i];  
  173. }
  174.       }
  175.     }
  176.     sample_f = sample_f/2;
  177.     actual_m = actual_m*2;
  178.     actual_n = actual_n*2;
  179.   }
  180. }
  181. #ifdef __STDC__
  182. bpconv(double *x_out, int lx, double *g0, double *g1, int lh,
  183.        double *x_inl, double *x_inh)
  184. #else
  185. bpconv(x_out, lx, g0, g1, lh, x_inl, x_inh)
  186. double *x_inl, *x_inh, *g0, *g1, *x_out;
  187. int lx, lh;
  188. #endif
  189. {
  190.   int i, j;
  191.   double x0;
  192.  
  193.   for (i=lh-2; i > -1; i--){
  194.     x_inl[i] = x_inl[lx+i];
  195.     x_inh[i] = x_inh[lx+i];
  196.   }
  197.   for (i=0; i<lx; i++){
  198.     x0 = 0;
  199.     for (j=0; j<lh; j++)
  200.       x0 = x0 + x_inl[j+i]*g0[lh-1-j] +
  201. x_inh[j+i]*g1[lh-1-j];
  202.     x_out[i] = x0;
  203.   }
  204. }