Trying to follow the math behind wavelets

Frnak McKenney wrote:

> So the "range" of the CWT(), the "surface" the CWT describes
> above the scale-time plane, is the same as the "range" of the
> signal(s), that is, amplitude.

That's correct.

> But the important thing (to me) is that I'm "pushing things
> around" and seeing results that match what my "model" predicts.

Right on!

> Hope y'all can enjoy it from afar.

It's indeed exciting (though something of an acquired taste), thanks 
for the explanation. See that everyone gets through the weather OK!


Martin

-- 
A trick that is used more than once is called a method.
--Ron Getoor paraphrasing George Polya and Gabor Szego

Martin Eisenberg wrote:

   ...

> [I] See that everyone gets through the weather OK!

I hope so! No problem here that I know of.

To guard against my canoe blowing around in the storm, I set it right 
side up and filled it with about half a ton of water. (I needn't have 
bothered.) The rain added five or six inches of additional water.

Jerry
-- 
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������

>
> Me, I'm still trying to understand the rules governing the selection
> of an appropriately-sized... �er, -scaled wavelet-screwdriver for
> pounding in a given set of signal-nails. &#4294967295;<grin!>
>
> Frank
> --

Hello Frank,

Have you tried Ingrid Daubechies' book "Ten Lectures on Wavelets."

She covers the rules behind wavelets, i.e., their being a complete
basis for example. Their construction ect.

Clay

Jerry Avins wrote:
> Martin Eisenberg wrote:

>> [I] See that everyone gets through the weather OK!

What I meant was not "I perceive" but "You take care". Should I have 
written it differently to convey that?


Martin

-- 
Seek simplicity and mistrust it.
--Alfred Whitehead

Martin Eisenberg wrote:
> Jerry Avins wrote:
>> Martin Eisenberg wrote:
> 
>>> [I] See that everyone gets through the weather OK!
> 
> What I meant was not "I perceive" but "You take care". Should I have 
> written it differently to convey that?

No. I jumped to a wrong conclusion. Good writers shape their sentences 
not merely to convey their intended meaning, but to make it impossible 
that any other meaning be attached. That's too high a standard for 
casual -- if heartfelt -- conversation. We need good readers, too.

Jerry
-- 
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������

On Sun, 7 Sep 2008 09:27:02 -0700 (PDT), clay@claysturner.com <clay@claysturner.com> wrote:

> Have you tried Ingrid Daubechies' book "Ten Lectures on Wavelets."
> 
> She covers the rules behind wavelets, i.e., their being a complete
> basis for example. Their construction ect.

Clay,

Sorry I didn't reply sooner.  I did see a copy of her book at the
University of Richmond Library, but I had trouble following more
than the basics of what she was doing.  Hubbard's "World of
Wavelets" is a bit more my speed, I'm afraid.  But thanks for the
suggestion.


Frank
--
    Little is inevitable in history, and even less in warfare. God
    may be thought to be on the side with the largest battalions,
    and quite often it works out in that manner, but wars often turn
    on an untold number of factors, including strategic and tactical
    decisions, good or bad leadership, the experience of the soldiery,
    the availability of resources, success in getting supplies where
    and when they are needed, the conduct of foreign policy,
    ideological zeal, willpower, and a measure of good luck.
                             -- John Ferling / Almost a Miracle
--
Frank McKenney, McKenney Associates
Richmond, Virginia / (804) 320-4887
Munged E-mail: frank uscore mckenney ayut mined spring dawt cahm (y'all)

On the off chance that someone muight be curious, I thought I'd post
a sample of a WWV waveform being "decoded" using wavelets.  It's a
simple pair of BCD digits modulating a 100Hz sine wave (assuming my
code is doing what I think it's doing <grin!>)

This is, of course, a nice, clean, theoretical signal; add a lot of
noise and it looks a lot messier.  That part is left as an E.F.S.
<grin!>

You'll need Scilab 4.1.2 and the Scilab Wavelet Toolbox (SWT) to run
the code.


Frank
------------------snip here-------------------------
//
// Scilab WWV and Wavelets Demo
//
// Requires Scilab Wavelet Toolbox
// Tested with WinXX Scilab v4.1.2, SWT v0.1.0-rc6.
//
// Written by: Frank McKenney, McKenney Associates
// Last updated: 08Oct2008
//
stacksize(floor(12*1000000));//      increase memory size
//
// Plotting Performance
// On my Dell 2.2GHz laptop, plotting more than about 20K points
// makes plotting slow (see also: non-interactive) and scaling and
// rotating the result equally sluggish.
// 
PlotMaxPoints = 30000;  // More than this makes plotting slow
//
// WindowSize = number of seconds of samples to fetch
// WindowOffset is offset into file in seconds
WindowSize = 10;                // Samples to fetch (seconds)
WindowOffset = 38;             // offset into file (seconds)
WindowPhaseShift = 0.30;  // Offset to put ticks on sec. boundary
//
// Add signals
//
Fs = 2500;  // "Sample" at 2500Hz
F0 = [100, 1.0]; // Hz, amp. Set to 0 to eliminate
F1 = [200, 1.0];   // Hz, amp. Set to 0 to eliminate
//
// Per NIST Time and Frequency Services 1394.pdf
//
// Time Signals: BCD code on 100 Hz subcarrier, 1 pulse/s
//
// "Bits are transmitted on the 100 Hz subcarrier using amplitude
//  modulation. A 200 ms pulse (20 cycles of 100 Hz) is used to
//  represent a 0 bit, and a 500 ms pulse (50 cycles of 100 Hz) is
//  used to represent a 1 bit. However, tone suppression deletes
//  the first 30 ms of each pulse. Therefore, 170 ms pulses are
//  recognized as 0 bits, and 470 ms pulses are recognized as 1 bits.
//  The leading edge of each pulse can serve as an on time marker,
//  but due to the tone suppression it actually occurs 30 ms after
//  the start of the second."
//
// Implications:
//
//  1) Fs needs to resolve 100Hz.
//  2) Fs needs to let us distinguish between 200/170ms trains of
//      100Hz and 500/470ms trains of 100Hz.
//
// Available Wavelet types, as of v0.1.0-rc6 (The Help file listing
// contains some name errors.)
//
//   haar  mexh  morl  cmor  cauchy  poisson
//   gaus1-8  cgau1-8  legd1-legd9
//   db1-db20  coif1-coif5  sym2-sym20
//   dmey  vaidyanathan  bior1.1-bior6.8
//   sinus  DOG  shan  fbsp
//   beylkin  bath4.0-bath4.15  bath6.0-bath6.15
//  
//  Scale choices will change with wavelet type
WaveletType = "beylkin"; // [db7 1:60], 
Scales = 1:60; // 1:12; // 1:30; 
//
// WWV sends BCD digits LSB first
BCDBits = [1 1 0 0  1 0 0 1]; // 93 in BCD, 8 bits = 8 seconds
//BCDBits = [1 0 1]; // 5 in BCD, 3 bits = 3 seconds
//BCDBits = [1]; // Single bit for testing
//BCDBits = [];   // No BCD signal
//
HalfSecond100Hz = sin( (F0(1)*2*%pi) * ([0:fix(Fs/2)-1]/Fs)   );
HalfSecond100Hz(1:fix( 0.030 * Fs ) ) = 0; // First 30msec silenced
Bits = HalfSecond100Hz;
HalfSecond100Hz( [1+fix( 0.200 * Fs ):fix(Fs/2)] ) = 0; 
Bits = [HalfSecond100Hz' , Bits' ]'; // Ugly, but laminates new row
Bits = [Bits, zeros(Bits)];
//
BCDSignal = Bits((1+BCDBits),:)'  ; // Select rows per index and ravel
BCDSignal = BCDSignal(:)'   ;      //  result (with transposes).
//
// Remove every other-than-nth value of s
//
function z=Decimate(s,n)
  z = s(1:n:$);
endfunction
//
// DEMO.  We're not really reading from a file here
//wf = "WWV-Richmond-2007-10-23-c2016-excellent-to-good-113sec-clip.wav";
//[WFcs, WFFs, WFbits] = wavread(wf,'size');  // Get WAV file information
//chs = WFcs(1); // Channel count
//ns = WFcs(2);   // Sample count
WFFs = Fs;
// Load samples
WSzSa = WindowSize * WFFs;                // Size (in samples)
// Offset in samples
WOSa = 0;
//WFy = wavread(wf,WSzSa+WOSa);
//WFy = WFy(WOSa+1:length(WFy));
//
DF = 1; // Reduce data size by Downsampling Factor
//y = Decimate(WFy,DF);   Fs = WFFs/DF;  // Decimate, adjust Fs
//WFy = []; //Free up storage
y = zeros(1,WindowSize * Fs); // Demo - no signal yet.
ny = length(y);
x = [1:ny];
//
// Add simulated WWV audio, matched to window
//
BCDsize = size(BCDSignal);  // No recursive indexing
if (BCDsize(1)) > 0 then;
  if (BCDsize(2) > length(y)) then;
    BCDSignal = BCDSignal(1:length(y)); // Too long - truncate
  elseif (BCDsize(2) < length(y)) then;
    BCDSignal = [BCDSignal, zeros(1,length(y)-length(BCDSignal))];
  end;
  y = y + BCDSignal;
end;
//
// Possibly amplify signal
//
//y = y * 3;
//
//  Possibly add reference signal(s) as markers
//
if (F0(1) > 0) then;
    yRef = F0(2) *  sin( (F0(1) *2*%pi) * ([0:length(y)-1]/Fs) );
    y = y + yRef;
end;
if (F1(1) > 0) then;
    yRef = F1(2) *  sin( (F1(1) *2*%pi) * ([0:length(y)-1]/Fs) );
    y = y + yRef;
end;
//
coef=cwt(y, Scales, WaveletType);
//
// "Downsample" the coefficients to a usable size array
//
[pts_y,pts_x] = size(coef);  // Note order here.
PlotDS = floor( (pts_x*pts_y) / PlotMaxPoints );
keepcols = 1:PlotDS:pts_x;
DScoef = coef(:,keepcols); // Reduce array size via columns (x)
//
// Plot CWT coefficients
//
set("figure_style","new");
clf();
f = gcf();
f.visible = "off"; // Eliminate frustrating redraws
f.figure_name = "Wavelet Analysis of WWV Signal";
set(f, "color_map", graycolormap(128));  //15 good
//
[xxx,c_cols] = size(DScoef);
x_time = WindowOffset + PlotDS * (0:c_cols-1)/Fs; // Time in seconds
//
DSc1 = DScoef  / max(DScoef);  // Normalize
DSc2 = DSc1 .* DSc1; DSc2 = DSc2 / max(DSc2);
mesh(x_time, Scales,  (DSc1 .* DSc1) );
//surf(x_time, Scales, (DScoef .* DScoef));
//plot3d1(x_time, Scales, (DScoef' .* DScoef'));
// Warning: grayplot() looks nice, but takes a long time
//grayplot(x_time, Scales, (DScoef .* DScoef)');
//
//
a=gca();
title = "Wavelet Analysis of WWV Signal";
title = title + sprintf("   Wavelet Type: %s",WaveletType);
title = title + sprintf("       Fs: %d   Ns: %d", Fs, WindowSize*Fs);
a.title.text = title;
a.x_label.text="Time (sec)";
a.y_label.text="Scale";
a.z_label.text="CWT()";
//
//a.rotation_angles = [6, -104]; // Viewpoint angles: alpha, theta
a.rotation_angles = [16, -108]; // Viewpoint angles: alpha, theta
//
// Draw contents of plot window
//
f.visible = "on";
//
//=============================================

--
    ...[T]he innovator has for enemies all those who have done well
    under the old conditions, and lukewarm defenders in those who
    may do well under the new.  -- Niccolo Machiavelli / The Prince
--
Frank McKenney, McKenney Associates
Richmond, Virginia / (804) 320-4887
Munged E-mail: frank uscore mckenney ayut mined spring dawt cahm (y'all)

Frnak McKenney wrote:

> On the off chance that someone muight be curious, I thought I'd
> post a sample of a WWV waveform being "decoded" using wavelets. 

So you haven't drowned after all. With hindsight, do you see any 
benefits of wavelets over other approaches to this task?


Martin

-- 
Quidquid latine scriptum est, altum videtur.

Martin,

Thanks for the interest.  I haven't given up on my research, but
other things keep interrupting.

On 11 Oct 2008 19:03:50 GMT, Martin Eisenberg <martin.eisenberg@udo.edu> wrote:
> Frnak McKenney wrote:
> 
>> On the off chance that someone muight be curious, I thought I'd
>> post a sample of a WWV waveform being "decoded" using wavelets. 
> 
> So you haven't drowned after all.

No, but we haven't had much rain lately, so there's still hope.
<grin!>

> ... With hindsight, do you see any 
> benefits of wavelets over other approaches to this task?

The "demo" I posted shows me that I _can_ use wavelets to clearly
pick WWV's digital time-code information out of an ideal (e.g.
Scilab-generated) signal.

On the other hand, using the same techniques on parts of a WAV file
recorded off a Heathkit MAC-II with an indoor antenna shows me that
this data is a lot harder to pick out from a RealWorld(tm) signal.
The audible portions of the recording are moderately clear most of
the time, so I'm assuming I didn't lose too much of the 100Hz during
recording, and I _think_ I can spot the data bits here and there,
but they come and go as I work my way through the WAV file.

Next step:  use wavelets to try reducing the noise.

Several of the papers in my collection describe what seems like a
reasonable technique:  dropping the signal gain at points where
there is "wideband" activity, that is, the wavelet "filtering" shows
presence simultaneously across a wide range of scales.

Hm.  And, just for grins (since it only now occurred to me), I think
I'll check the 100Hz performance of the recording system's audio
support; all the K8N-E book says is "Realtek ALC850", so I may have
to dig out my signal generator.


Frank McKenney, McKenney Associates
Richmond, Virginia / (804) 320-4887
Munged E-mail: frank uscore mckenney ayut mined spring dawt cahm (y'all)