zhangyiss/wavelet1d
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Periodic DWT Code Modified

Browse Source 
ChangeLog ver 0.3

A couple of major changes and a few minor tweaks.

1. Periodic DWT code changed. The issue was Perfect Reconstruction wasn't working when filter length became larger than
 decmated signal. Code changed to implement actual periodic extension instead of using circular shift to mimic periodic 
extension. Change is strictly under the hood and header file remains unchanged.

2. Dyadic Zero padding is changed to padding the signal by the last value of the signal instead of zeros.
This commit is contained in:
rafat.hsn@gmail.com 2011-05-28 03:20:59 +00:00
parent d26770e025
commit efc1d49cfe
1 changed files with 5258 additions and 5302 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

206
src/wavelet.cpp
View File

@ -441,7 +441,7 @@ void* iswt(vector<double> &swtop,int J, string nm, vector<double> &iswt_output)
void* swt(vector<double> &sig, int J, string nm, vector<double> &swt_output, int &length) { void* swt(vector<double> &sig, int J, string nm, vector<double> &swt_output, int &length) {
vector<double> lpd, hpd, lpr, hpr; vector<double> lpd, hpd, lpr, hpr;
dyadic_zpad_1d(sig); // dyadic_zpad_1d(sig);
int N = sig.size(); int N = sig.size();
length = N; length = N;
@ -1738,18 +1738,19 @@ void* dwt(vector<double> &sig, int J, string nm, vector<double> &dwt_output
, vector<double> &flag ) { , vector<double> &flag ) {
int Max_Iter; int Max_Iter;
Max_Iter = (int) ceil(log( double(sig.size()))/log (2.0)); Max_Iter = (int) ceil(log( double(sig.size()))/log (2.0)) - 1;
if ( Max_Iter < J) { if ( Max_Iter < J) {
cout << J << " Iterations are not possible with signal of length " << sig.size() << endl; cout << J << " Iterations are not possible with signal of length " << sig.size() << endl;
exit(1); exit(1);
} }
vector<double> original_copy, appx_sig, det_sig; vector<double> original_copy,orig, appx_sig, det_sig;
original_copy = sig; original_copy = sig;
// Zero Pad the Signal to nearest 2^ M value ,where M is an integer. // Zero Pad the Signal to nearest 2^ M value ,where M is an integer.
unsigned int n_size = sig.size(); unsigned int n_size = sig.size();
flag.push_back(0); flag.push_back(0);
double int_n_size = log10 (static_cast<double> (n_size)) / log10(2.0); double int_n_size = log10 (static_cast<double> (n_size)) / log10(2.0);
if ( (pow(2.0, double(int(ceil(int_n_size)))) - pow(2.0, int_n_size)) != 0) { if ( (pow(2.0, double(int(ceil(int_n_size)))) - pow(2.0, int_n_size)) != 0) {
dyadic_zpad_1d(sig); dyadic_zpad_1d(sig);
@ -1759,6 +1760,7 @@ void* dwt(vector<double> &sig, int J, string nm, vector<double> &dwt_output
} }
flag.push_back(J); flag.push_back(J);
orig = sig;
ofstream appx("appx.txt", ios::trunc); ofstream appx("appx.txt", ios::trunc);
appx.close(); appx.close();
@ -1769,8 +1771,9 @@ void* dwt(vector<double> &sig, int J, string nm, vector<double> &dwt_output
// Storing Signal for Gnuplot // Storing Signal for Gnuplot
unsigned int count = sig.size(); unsigned int count = sig.size();
ofstream gnusig("gnusig.dat"); ofstream gnusig("gnusig.dat");
for (unsigned int i = 0;i < count; i++) { unsigned int ct_orig=original_copy.size();
gnusig << i << " " << sig[i] << endl; for (unsigned int i = 0;i < ct_orig; i++) {
gnusig << i << " " << original_copy[i] << endl;
} }
gnusig.close(); gnusig.close();
@ -1787,7 +1790,7 @@ void* dwt(vector<double> &sig, int J, string nm, vector<double> &dwt_output
for (int iter = 0; iter < J; iter++) { for (int iter = 0; iter < J; iter++) {
dwt1(nm,sig, appx_sig, det_sig); dwt1(nm,orig, appx_sig, det_sig);
dwt_output.insert(dwt_output.begin(),det_sig.begin(),det_sig.end()); dwt_output.insert(dwt_output.begin(),det_sig.begin(),det_sig.end());
@ -1799,8 +1802,8 @@ void* dwt(vector<double> &sig, int J, string nm, vector<double> &dwt_output
// vector<double> lp1,hp1,lp2,hp2; // vector<double> lp1,hp1,lp2,hp2;
// filtcoef(nm,lp1,hp1,lp2,hp2); // filtcoef(nm,lp1,hp1,lp2,hp2);
// gnudwtplot(sig, appx_sig, det_sig,lp1,hp1,lp2,hp2); // gnudwtplot(sig, appx_sig, det_sig,lp1,hp1,lp2,hp2);
sig.clear(); // sig.clear();
sig = appx_sig; orig = appx_sig;
appx_sig.clear(); appx_sig.clear();
det_sig.clear(); det_sig.clear();
@ -1879,7 +1882,7 @@ void* dwt(vector<double> &sig, int J, string nm, vector<double> &dwt_output
} }
void circshift(vector<double> &sig_cir, int L){ void circshift(vector<double> &sig_cir, int L){
if ( abs(L) > sig_cir.size()) { if ( abs(L) >(signed int) sig_cir.size()) {
L = sign(L) * (abs(L) % sig_cir.size()); L = sign(L) * (abs(L) % sig_cir.size());
} }
@ -1961,8 +1964,9 @@ void* dwt1(string wname, vector<double> &signal, vector<double> &cA, vector<doub
int len_lpfilt = lpd.size(); int len_lpfilt = lpd.size();
int len_hpfilt = hpd.size(); int len_hpfilt = hpd.size();
int len_avg = (len_lpfilt + len_hpfilt) / 2; int len_avg = (len_lpfilt + len_hpfilt) / 2;
// cout << len_lpfilt << "Filter" << endl; // cout << len_lpfilt << "Filter" << endl;
circshift(signal,-len_avg / 2); // Signal is shifted circularly in order to perform per_ext(signal,len_avg / 2); // Periodic Extension
// computations designed to deal with boundary distortions // computations designed to deal with boundary distortions
// Low Pass Filtering Operations in the Analysis Filter Bank Section // Low Pass Filtering Operations in the Analysis Filter Bank Section
@ -1974,14 +1978,10 @@ void* dwt1(string wname, vector<double> &signal, vector<double> &cA, vector<doub
// Downsampling by 2 gives cA // Downsampling by 2 gives cA
downsamp(cA_undec, D, cA); downsamp(cA_undec, D, cA);
int L = len_lpfilt / 2;
int N = len_sig / 2; int N = len_sig / 2;
for (int i = 0; i < L; i++) { cA.erase(cA.begin(),cA.begin()+len_avg/2);
cA[i] = cA[N+i]+cA[i]; cA.erase(cA.end()-len_avg/2,cA.end());
}
cA.resize(N);
circshift(cA,len_avg/2 );
// High Pass Filtering Operations in the Analysis Filter Bank Section // High Pass Filtering Operations in the Analysis Filter Bank Section
// int len_cA =(int) floor(double (len_sig + len_lpfilt -1) / double (2)); // int len_cA =(int) floor(double (len_sig + len_lpfilt -1) / double (2));
@ -1994,13 +1994,8 @@ void* dwt1(string wname, vector<double> &signal, vector<double> &cA, vector<doub
// Downsampling by 2 gives cA // Downsampling by 2 gives cA
downsamp(cD_undec, D, cD); downsamp(cD_undec, D, cD);
int Lh = len_hpfilt / 2; cD.erase(cD.begin(),cD.begin()+len_avg/2);
cD.erase(cD.end()-len_avg/2,cD.end());
for (int i = 0; i < Lh; i++) {
cD[i] = cD[N+i]+cD[i];
}
cD.resize(N);
circshift(cD,len_avg/2 );
// File Outputs // File Outputs
@ -2024,11 +2019,12 @@ void* dyadic_zpad_1d(vector<double> &signal) {
int D = (int) ceil(M); int D = (int) ceil(M);
double int_val = pow(2.0, double(D)) - pow(2.0, M); double int_val = pow(2.0, double(D)) - pow(2.0, M);
int zeros = (int) int_val; int z = (int) int_val;
vector<double>::iterator a_it; vector<double>::iterator a_it;
a_it = signal.end(); a_it = signal.end();
// double val = signal[N-1]; double val = signal[N-1];
signal.insert(a_it,zeros,0); // double val = 0;
signal.insert(a_it,z,val);
return 0; return 0;
} }
@ -2087,9 +2083,8 @@ void* idwt(vector<double> &dwtop,vector<double> &flag, string nm,
vector<double> app; vector<double> app;
vector<double> detail; vector<double> detail;
unsigned int app_len = dwtop.size() / int(pow(2.0,J)); unsigned int app_len = dwtop.size() / int(pow(2.0,J));
vector<double>::iterator dwt;
dwt = dwtop.begin(); app.assign(dwtop.begin(),dwtop.begin()+app_len);
app.assign(dwt,dwtop.begin()+app_len);
detail.assign(dwtop.begin()+app_len, dwtop.begin()+ 2* app_len); detail.assign(dwtop.begin()+app_len, dwtop.begin()+ 2* app_len);
for (int i = 0; i < J; i++) { for (int i = 0; i < J; i++) {
@ -2133,14 +2128,14 @@ void* idwt1(string wname, vector<double> &X, vector<double> &cA, vector<double>
int len_lpfilt = lpr1.size(); int len_lpfilt = lpr1.size();
int len_hpfilt = hpr1.size(); int len_hpfilt = hpr1.size();
int len_avg = (len_lpfilt + len_hpfilt) / 2; int len_avg = (len_lpfilt + len_hpfilt) / 2;
unsigned int N = 2 * cA.size(); //unsigned int N = 2 * cA.size();
int U = 2; // Upsampling Factor int U = 2; // Upsampling Factor
// Operations in the Low Frequency branch of the Synthesis Filter Bank // Operations in the Low Frequency branch of the Synthesis Filter Bank
vector<double> cA_up; vector<double> cA_up;
vector<double> X_lp; vector<double> X_lp;
circshift(cA,-len_avg/2); per_ext(cA,len_avg/2);
upsamp(cA, U, cA_up); upsamp(cA, U, cA_up);
convfft(cA_up, lpr1, X_lp); convfft(cA_up, lpr1, X_lp);
@ -2150,39 +2145,21 @@ void* idwt1(string wname, vector<double> &X, vector<double> &cA, vector<double>
vector<double> cD_up; vector<double> cD_up;
vector<double> X_hp; vector<double> X_hp;
circshift(cD,-len_avg/2); per_ext(cD,len_avg/2);
upsamp(cD, U, cD_up); upsamp(cD, U, cD_up);
convfft(cD_up, hpr1, X_hp); convfft(cD_up, hpr1, X_hp);
// Operations to obtain reconstructed signal by folding back and circular shifting
for (unsigned int i = 0 ; i < X_lp.size() - N ; i++){
X_lp[i] = X_lp[N+i] + X_lp[i];
}
for (unsigned int i = 0 ; i < X_hp.size() - N; i++){
X_hp[i] = X_hp[N+i] + X_hp[i];
}
X_lp.resize(N);
X_hp.resize(N);
// for (unsigned int i =0 ; i < N ; i++){
// X[i] = X_lp[i] + X_hp[i];
// }
vecsum(X_lp,X_hp,X); vecsum(X_lp,X_hp,X);
// Work on circular shift // Remove periodic extension
circshift(X,len_avg +len_avg / 2 -1);
X.erase(X.begin(),X.begin()+len_avg+len_avg/2-1);
X.erase(X.end()-len_avg-len_avg/2,X.end());
/*
for (unsigned int i = 0; i < N; i++){
sig << X[i] << endl;
}
*/
return 0; return 0;
} }
@ -2521,96 +2498,76 @@ void* branch_lp_hp_up(string wname,vector<double> &cA, vector<double> &cD, vecto
int len_lpfilt = lpr1.size(); int len_lpfilt = lpr1.size();
int len_hpfilt = hpr1.size(); int len_hpfilt = hpr1.size();
int len_avg = (len_lpfilt + len_hpfilt) / 2; int len_avg = (len_lpfilt + len_hpfilt) / 2;
unsigned int N = 2 * cA.size(); //unsigned int N = 2 * cA.size();
int U = 2; // Upsampling Factor int U = 2; // Upsampling Factor
// Operations in the Low Frequency branch of the Synthesis Filter Bank // Operations in the Low Frequency branch of the Synthesis Filter Bank
vector<double> cA_up; vector<double> cA_up;
vector<double> X_lp; vector<double> X_lp;
circshift(cA,-len_avg/2); per_ext(cA,len_avg/2);
upsamp(cA, U, cA_up); upsamp(cA, U, cA_up);
convfft(cA_up, lpr1, X_lp); convfft(cA_up, lpr1, X_lp);
// Operations in the High Frequency branch of the Synthesis Filter Bank // Operations in the High Frequency branch of the Synthesis Filter Bank
vector<double> cD_up; vector<double> cD_up;
vector<double> X_hp; vector<double> X_hp;
circshift(cD,-len_avg/2); per_ext(cD,len_avg/2);
upsamp(cD, U, cD_up); upsamp(cD, U, cD_up);
convfft(cD_up, hpr1, X_hp); convfft(cD_up, hpr1, X_hp);
// Operations to obtain reconstructed signal by folding back and circular shifting
for (unsigned int i = 0 ; i < X_lp.size() - N ; i++){
X_lp[i] = X_lp[N+i] + X_lp[i];
}
for (unsigned int i = 0 ; i < X_hp.size() - N; i++){
X_hp[i] = X_hp[N+i] + X_hp[i];
}
X_lp.resize(N);
X_hp.resize(N);
// for (unsigned int i =0 ; i < N ; i++){
// X[i] = X_lp[i] + X_hp[i];
// }
vecsum(X_lp,X_hp,X); vecsum(X_lp,X_hp,X);
// Work on circular shift // Remove periodic extension
circshift(X,len_avg +len_avg / 2 -1);
return 0; X.erase(X.begin(),X.begin()+len_avg+len_avg/2-1);
X.erase(X.end()-len_avg-len_avg/2,X.end());
return 0;
} }
void* branch_hp_dn(string wname, vector<double> &signal, vector<double> &sigop) { void* branch_hp_dn(string wname, vector<double> &signal, vector<double> &sigop) {
vector<double> lpd1,hpd1,lpr1,hpr1; vector<double> lpd, hpd, lpr, hpr;
filtcoef(wname,lpd1,hpd1,lpr1,hpr1);
filtcoef(wname,lpd,hpd,lpr,hpr);
// for (unsigned int i = 0; i < signal.size(); i++) {
// cout << signal[i] << endl;
// out2 << signal[i] <<endl;
// }
unsigned int temp_len = signal.size(); unsigned int temp_len = signal.size();
if ( (temp_len % 2) != 0) { if ( (temp_len % 2) != 0) {
double temp = 0; double temp =signal[temp_len - 1];
signal.push_back(temp); signal.push_back(temp);
} }
int len_sig = signal.size(); int len_lpfilt = lpd.size();
int len_hpfilt = hpd1.size(); int len_hpfilt = hpd.size();
int len_avg = len_hpfilt; int len_avg = (len_lpfilt + len_hpfilt) / 2;
// cout << len_lpfilt << "Filter" << endl; // cout << len_lpfilt << "Filter" << endl;
circshift(signal,-len_avg / 2); // Signal is shifted circularly in order to perform per_ext(signal,len_avg / 2); // Periodic Extension
// computations designed to deal with boundary distortions // computations designed to deal with boundary distortions
// Low Pass Filtering Operations in the Analysis Filter Bank Section
// int len_cA =(int) floor(double (len_sig + len_lpfilt -1) / double (2));
// High Pass Filtering Operations in the Analysis Filter Bank Section vector<double> cA_undec;
// int len_cA =(int) floor(double (len_sig + len_lpfilt -1) / double (2));
vector<double> sig_undec;
// convolving signal with lpd, Low Pass Filter, and O/P is stored in cA_undec // convolving signal with lpd, Low Pass Filter, and O/P is stored in cA_undec
convfft(signal,hpd1,sig_undec); convfft(signal,hpd,cA_undec);
// Downsampling Factor is 2 int D = 2; // Downsampling Factor is 2
int D =2 ;
// Downsampling by 2 gives cA // Downsampling by 2 gives cA
downsamp(sig_undec, D, sigop); downsamp(cA_undec, D, sigop);
int Lh = len_hpfilt / 2;
int N = len_sig / 2;
for (int i = 0; i < Lh; i++) {
sigop[i] = sigop[N+i]+sigop[i];
}
sigop.resize(N);
circshift(sigop,len_avg/2 );
filtcoef(wname,lpd1,hpd1,lpr1,hpr1);
sigop.erase(sigop.begin(),sigop.begin()+len_avg/2);
sigop.erase(sigop.end()-len_avg/2,sigop.end());
return 0; return 0;
@ -2618,43 +2575,42 @@ return 0;
} }
void* branch_lp_dn(string wname, vector<double> &signal, vector<double> &sigop){ void* branch_lp_dn(string wname, vector<double> &signal, vector<double> &sigop){
vector<double> lpd1,hpd1,lpr1,hpr1; vector<double> lpd, hpd, lpr, hpr;
filtcoef(wname,lpd1,hpd1,lpr1,hpr1); filtcoef(wname,lpd,hpd,lpr,hpr);
// for (unsigned int i = 0; i < signal.size(); i++) {
// cout << signal[i] << endl;
// out2 << signal[i] <<endl;
// }
unsigned int temp_len = signal.size(); unsigned int temp_len = signal.size();
if ( (temp_len % 2) != 0) { if ( (temp_len % 2) != 0) {
double temp = 0; double temp =signal[temp_len - 1];
signal.push_back(temp); signal.push_back(temp);
} }
int len_sig = signal.size(); int len_lpfilt = lpd.size();
int len_lpfilt = lpd1.size(); int len_hpfilt = hpd.size();
int len_avg = len_lpfilt; int len_avg = (len_lpfilt + len_hpfilt) / 2;
// cout << len_lpfilt << "Filter" << endl; // cout << len_lpfilt << "Filter" << endl;
circshift(signal,-len_avg / 2); // Signal is shifted circularly in order to perform per_ext(signal,len_avg / 2); // Periodic Extension
// computations designed to deal with boundary distortions // computations designed to deal with boundary distortions
// Low Pass Filtering Operations in the Analysis Filter Bank Section // Low Pass Filtering Operations in the Analysis Filter Bank Section
// int len_cA =(int) floor(double (len_sig + len_lpfilt -1) / double (2)); // int len_cA =(int) floor(double (len_sig + len_lpfilt -1) / double (2));
vector<double> sig_undec; vector<double> cA_undec;
// convolving signal with lpd, Low Pass Filter, and O/P is stored in cA_undec // convolving signal with lpd, Low Pass Filter, and O/P is stored in cA_undec
convfft(signal,lpd1,sig_undec); convfft(signal,lpd,cA_undec);
int D = 2; // Downsampling Factor is 2 int D = 2; // Downsampling Factor is 2
// Downsampling by 2 gives cA // Downsampling by 2 gives cA
downsamp(sig_undec, D, sigop); downsamp(cA_undec, D, sigop);
int L = len_lpfilt / 2;
int N = len_sig / 2;
for (int i = 0; i < L; i++) { sigop.erase(sigop.begin(),sigop.begin()+len_avg/2);
sigop[i] = sigop[N+i]+ sigop[i]; sigop.erase(sigop.end()-len_avg/2,sigop.end());
}
sigop.resize(N);
circshift(sigop,len_avg/2 );
filtcoef(wname,lpd1,hpd1,lpr1,hpr1);
return 0; return 0;