-
Notifications
You must be signed in to change notification settings - Fork 2
/
mdwt_r.c
executable file
·178 lines (159 loc) · 5.69 KB
/
mdwt_r.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
/*
File Name: MDWT.c
Last Modification Date: 06/14/95 13:15:44
Current Version: MDWT.c 2.4
File Creation Date: Wed Oct 19 10:51:58 1994
Author: Markus Lang <lang@jazz.rice.edu>
Copyright (c) 2000 RICE UNIVERSITY. All rights reserved.
Created by Markus Lang, Department of ECE, Rice University.
This software is distributed and licensed to you on a non-exclusive
basis, free-of-charge. Redistribution and use in source and binary forms,
with or without modification, are permitted provided that the following
conditions are met:
1. Redistribution of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
2. Redistribution in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
3. All advertising materials mentioning features or use of this software
must display the following acknowledgment: This product includes
software developed by Rice University, Houston, Texas and its contributors.
4. Neither the name of the University nor the names of its contributors
may be used to endorse or promote products derived from this software
without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY WILLIAM MARSH RICE UNIVERSITY, HOUSTON, TEXAS,
AND CONTRIBUTORS AS IS AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL RICE UNIVERSITY
OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
OR BUSINESS INTERRUPTIONS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
OTHERWISE), PRODUCT LIABILITY, OR OTHERWISE ARISING IN ANY WAY OUT OF THE
USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
For information on commercial licenses, contact Rice University's Office of
Technology Transfer at techtran@rice.edu or (713) 348-6173
Change History: Fixed the code such that 1D vectors passed to it can be in
either passed as a row or column vector. Also took care of
the code such that it will compile with both under standard
C compilers as well as for ANSI C compilers
Jan Erik Odegard <odegard@ece.rice.edu> Wed Jun 14 1995
%y = mdwt(x,h,L);
%
% function computes the discrete wavelet transform y for a 1D or 2D input
% signal x.
%
% Input:
% x : finite length 1D or 2D signal (implicitely periodized)
% h : scaling filter
% L : number of levels. in case of a 1D signal length(x) must be
% divisible by 2^L; in case of a 2D signal the row and the
% column dimension must be divisible by 2^L.
%
% see also: midwt, mrdwt, mirdwt
*/
#include <math.h>
#include <stdio.h>
#define max(A,B) (A > B ? A : B)
#define mat(a, i, j) (*(a + (m*(j)+i))) /* macro for matrix indices */
#ifdef __STDC__
MDWT(double *x, int m, int n, double *h, int lh, int L, double *y)
#else
MDWT(x, m, n, h, lh, L, y)
double *x, *h, *y;
int m, n, lh, L;
#endif
{
double *h0, *h1, *ydummyl, *ydummyh, *xdummy;
long i, j;
int actual_L, actual_m, actual_n, r_o_a, c_o_a, ir, ic, lhm1;
xdummy = (double *)mxCalloc(max(m,n)+lh-1,sizeof(double));
ydummyl = (double *)mxCalloc(max(m,n),sizeof(double));
ydummyh = (double *)mxCalloc(max(m,n),sizeof(double));
h0 = (double *)mxCalloc(lh,sizeof(double));
h1 = (double *)mxCalloc(lh,sizeof(double));
/* analysis lowpass and highpass */
if (n==1){
n = m;
m = 1;
}
for (i=0; i<lh; i++){
h0[i] = h[lh-i-1];
h1[i] =h[i];
}
for (i=0; i<lh; i+=2)
h1[i] = -h1[i];
lhm1 = lh - 1;
actual_m = 2*m;
actual_n = 2*n;
/* main loop */
for (actual_L=1; actual_L <= L; actual_L++){
if (m==1)
actual_m = 1;
else{
actual_m = actual_m/2;
r_o_a = actual_m/2;
}
actual_n = actual_n/2;
c_o_a = actual_n/2;
/* go by rows */
for (ir=0; ir<actual_m; ir++){ /* loop over rows */
/* store in dummy variable */
for (i=0; i<actual_n; i++)
if (actual_L==1)
xdummy[i] = mat(x, ir, i);
else
xdummy[i] = mat(y, ir, i);
/* perform filtering lowpass and highpass*/
fpsconv(xdummy, actual_n, h0, h1, lhm1, ydummyl, ydummyh);
/* restore dummy variables in matrices */
ic = c_o_a;
for (i=0; i<c_o_a; i++){
mat(y, ir, i) = ydummyl[i];
mat(y, ir, ic++) = ydummyh[i];
}
}
/* go by columns in case of a 2D signal*/
if (m>1){
for (ic=0; ic<actual_n; ic++){ /* loop over column */
/* store in dummy variables */
for (i=0; i<actual_m; i++)
xdummy[i] = mat(y, i, ic);
/* perform filtering lowpass and highpass*/
fpsconv(xdummy, actual_m, h0, h1, lhm1, ydummyl, ydummyh);
/* restore dummy variables in matrix */
ir = r_o_a;
for (i=0; i<r_o_a; i++){
mat(y, i, ic) = ydummyl[i];
mat(y, ir++, ic) = ydummyh[i];
}
}
}
}
}
#ifdef __STDC__
fpsconv(double *x_in, int lx, double *h0, double *h1, int lhm1,
double *x_outl, double *x_outh)
#else
fpsconv(x_in, lx, h0, h1, lhm1, x_outl, x_outh)
double *x_in, *h0, *h1, *x_outl, *x_outh;
int lx, lhm1;
#endif
{
int i, j, ind;
double x0, x1;
for (i=lx; i < lx+lhm1; i++)
x_in[i] = *(x_in+(i-lx));
ind = 0;
for (i=0; i<(lx); i+=2){
x0 = 0;
x1 = 0;
for (j=0; j<=lhm1; j++){
x0 = x0 + x_in[i+j]*h0[lhm1-j];
x1 = x1 + x_in[i+j]*h1[lhm1-j];
}
x_outl[ind] = x0;
x_outh[ind++] = x1;
}
}