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)

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.

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)

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)

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.
�����������������������������������������������������������������������

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

>
> 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

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.
�����������������������������������������������������������������������

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

Hi, Martin.

If I seem a bit more confused than usual it's because I'm
alternating writing my reply with a periodic background sampling
process:  "The Effects of Extended Hurricane-Based Rainfall on
Mostly-Impermeable Cinder Block Basement Walls".  The good news is
that my sampling device (a.k.a.  "ShopVac") has a 16 gallon buffer.

On 5 Sep 2008 15:33:20 GMT, Martin Eisenberg <martin.eisenberg@udo.edu> wrote:
> Frnak McKenney wrote:
> 
>>>> Assumption:  The CWT(g(),psi(),t,s) is, in some measure,
>>>> related to the Fourier "frequency" spectrum of f().  That is,
>>>> for a fixed g() and psi(), if the frequency spectrum of g()
>>>> contains some frequency f0 at time t0, then for some 'scale'
>>>> s0 the results of the CWT(g(),psi(),t0,s0) has some sort of
>>>> "peak".
>>> 
>>> It has power there but that doesn't say anything about the
>>> waveform carrying it.
>> 
>> "Power".  As in "magnitude squared".  Okay...  the signal g()
>> has amplitude-hence-magnitude, and a wavelet has amplitude-hence
>> -magnitude, so even though my head hurts when I try to visualize
>> it, I can accept the idea that g()*psi() -- or, more precisely,
>> the sum/integral of g()*psi() -- represents "power".
> 
> Argh no, erase that thought ;) It's nothing to do with the mechanics 
> of computation! The signal has some spectrotemporal distribution of 
> power as a matter of fact, and the scalogram is designed to present 
> that in a recognizable way (one way among others, of course).

Okay.  I misunderstood.  I thought I had a handle on what the range
(output) of a CWT(scale,time) was, and my initial reading of your
reply threw me a bit.

Let me try to put my current "model" into words and see if it makes 
sense:

The CWT(f(),phi(),s,tau) maps a single function/signal f() into a
_family_ of function/signals.  Each member of this (possibly
infinite) family is created by repeatedly taking, at each point of
f(), an inner product of f() with a scaled copy of the wavelet
function psi() which has been translated to that point.

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.
 
>> Another question, or at least request for confirmation.  Suppose
>> I have the (infinite) CWT scale-time (y-x) plane laid out in
>> front of me, and further assume that I can stretch out a pair of
>> more-than- infinitely-extensible hands to squeeze it from the
>> top and bottom (Ack!)  into a horizontal line.  If my chosen
>> wavelets form an orthonormal basis (that is, they chop up
>> functions in such a way that they can be exactly reconstructed),
>> don't I get back my original function?
> 
> You can imagine it that way, noting that the synthesis wavelet 
> (reconstruction filterbank) gets thrown in too.

Oh.  Right.  I forgot that the "chopping up" involved changing the
shape of the signal "pieces" in a specific fashion (dependent on the
choice of wavelet), and therefore one would have to "undo" those
changes in order to put the "pieces" back together properly.

Hey!  That actually made sense!  (And who said miracles were out of
fashion?  <grin!>
 
>> Okay, that seems fairly trivial.  But...  suppose I squeeze from
>> the right and left?  Well, if I've left in any unbounded
>> thingies, such as sine-waves-extended-to-forever, some or all of
>> the result may blow up in my face. But ignoring that messy (but
>> very real) possibility, it "feels" like I wind up with a set of
>> values representing the "scales" at which my original function
>> had "power" at any point in time.
>> 
>> What would you call it?  A...  "scale spectrum"?
> 
> Good question; I don't know if there's an established term. But you 
> could also describe it without reference to wavelets as a Fourier 
> spectrum aggregated over log-spaced intervals.

After listening to your restatement, I think I just re-invented the
Constant-Q Transform in "scale" clothing. <grin!>

But the important thing (to me) is that I'm "pushing things around"
and seeing results that match what my "model" predicts.
 
So now I go back and re-read Hubbard and the five or six papers that
seemed to make the most sense.  If I'm lucky, a few more bits will
make sense this time around.
 
>> Meanwhile, in TheRealWorld(tm), it appears that Murphy -- long
>> recognized as the patron saint of Data Processing -- has decided
>> to take a hand in the American Presidential Election process. 
>> "May we live in interesting times", indeed.  <grin!>
> 
> I'm not American, you know. Does that mean you consider the election 
> decided with McCain's nomination?

Not in the least. But let me explain, or at least try to. <grin!>

Every four years we hold our Presidential elections.  As a part of
the process each candidate must undergo a long series of trials to
show that they are Worthy:  Trial by Insult, Trial by Shaming, Trial
by Prolonged Examination in Minute Detail (a.k.a.  Trial by Media),
Trial of Ideological Purity, Trial of Historical Consistency, Trial
by Slings and Arrows of Outrageous Fortune, ...  well, you get the
idea.  Points are deducted each time a candidate fails to stoically
endure one of their trials; points are added when they can not
merely endure a trial better/longer than an opponent, but make it
appear as if the results are unimportant to them.
 
Until a week ago (only one week?)  it appeared that we would be
witness to the trials of four fairly similar candidates.  The three
"knowns" were all legislators (two old-timers and one newcomer), and
most of us assumed that McCain's choice for VP would be yet another
Washington "insider".

McCain's selection of Governor Sarah Palin (Alaska) as his running
mate has put new life into the process.  She's not well known, so
the media can spend hours discussing her without getting bored.
She's female, as governor she has experience none of the other three
have, she's conservative, she hunts and fishes, and she has
successfully fought the oil companies, so the political pundits have
plenty of material for speculating on (e.g.)  what percentage of
Senator Clinton's former supporters will cross over to vote for her.
She's attractive and personable enough, and can speak coherent
sentences, so she's likely to do well on talk shows.  Finally, her
personality and views have enough "edges" for easy caricature, which
makes her a prime candidate for (e.g.)  Saturday Night Live.
<grin!>

Even the people who don't like her are talking about her.

As I say, "interesting times".  What looked to be a fairly
predictable (hence boring) see-saw battle between the parties for
the last few independent voters has suddenly shifted into unexpected
dimensions, and this is likely to increase the voter turnout on both
sides.  And since that puts the outcome of the election even more in
doubt, people (voters) are paying even more attention to the
process.  <grin!>

Anyway, that's my not-so-humble opinion. Hope y'all can enjoy it
from afar.


Frank
--
    The Democrats are the party of government activism, the party
    that says government can make you richer, smarter, taller, and
    get the chickweed out of your lawn. Republicans are the party
    that says government doesn't work, and then get elected and
    prove it.                -- P. J. O'Rourke
--
Frank McKenney, McKenney Associates
Richmond, Virginia / (804) 320-4887
Munged E-mail: frank uscore mckenney ayut mined spring dawt cahm (y'all)