-
Notifications
You must be signed in to change notification settings - Fork 0
/
HaarDWT.c
91 lines (62 loc) · 1.52 KB
/
HaarDWT.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
#include "HaarDWT.h"
#include <stdio.h>
#include <string.h>
const float sqrthalf = 0.70710678118;
static void rec_coef_1D_HaarDWT(float *ip, float *op, long n);
void rec_coef_1D_HaarDWT(float *ip, float *op, long n){
if(n==1)
return; // Base case
int i = 0;
for ( ; i<n/2; ++i) // approx, scale, lo freq
op[i] = ip[2*i] + ip[2*i+1];
for ( ; i<n; ++i) // detail, wavelet, hi freq
op[i] = ip[2*i-n] - ip[2*i-n+1];
for ( i=0 ; i<n; ++i) // normailze
op[i] *= sqrthalf;
// for ( i=0 ; i<n; ++i)
// printf("%f\t", op[i]);
// printf("\n");
float ip2[n/2];
for ( i=0 ; i<n/2; ++i)
ip2[i] = op[i];
rec_coef_1D_HaarDWT(ip2, op, n/2);
}
void coef_1D_HaarDWT(float *ip, float *op, long n){
if( (n==0) || (n&(n-1)) ){ // Check if n is not a power of 2
printf("%ld is not a power of 2\n", n);
return;
}
rec_coef_1D_HaarDWT(ip, op, n);
}
void coef_1D_if_HaarDWT(int *ip, float *op, long n){
if( (n==0) || (n&(n-1)) ){ // Check if n is not a power of 2
printf("%ld not a power of 2\n", n);
return;
}
if(n==1){
op[0] = ip[0];
return;
}
int tp[n];
memcpy(tp, ip, sizeof(tp)); // Copy array ip to tp
for(int k=1; k<n; k*=2){ // k=2^j, j is scale
for(int i=0; i<n; i+=2*k){
tp[i] = tp[i] + tp[i+k];
tp[i+k] = tp[i] - 2*tp[i+k];
}
}
int j = n-1;
for(int k=1; k<n; k*=2){
for(int i=k; i<n; i+=2*k, j--){
op[j] = tp[i];
}
}
op[0] = tp[0];
for(int k=2; k<=n; k*=2){
for(int i=0; i<k; ++i)
op[i] *= sqrthalf;
}
}
void sig_1D_HaarDWT(float *ip, float *op, long n){
return;
}