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)
WindowPhase****ft = 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)
|
