an initial experiment, Aug 4 2023

Transmitter RWM in Moscow transmits a CW waveform for around 8 minutes on the half hour, at frequencies that are close to WWV, namely (see https://en.wikipedia.org/wiki/RWM)  4.996 MHz with 5 kW and on 9.996 and 14.996 MHz with 8 kW.^[1]

^{Use a web-SDR receiver in the UK to capture this signal.  It's http://thescotts.ddns.net:8073/.} 

^{data file: thescotts.ddns.net_2023-08-04T05_59_28Z_9996.08_cwn.wav}

^{data were captured at 12KHz and the detuning gave a beat frequency of around 500 Hz so well into the Nyquist regime. 
That gives us insight into propagation variations at the 200 Hz or so level, around 5 mseec. In general if we have an SDR with some instantaneous BW Df, we can detune to get the beat at DR/2 and so 2/DF is the temporal resolution of looking at phase variations.} 

^{beat frequency is around 400 Hz, beat period is therefore 2.5 msec and there are around 24 points per beat cycle.} 

^{I've convinced myself that a phase change of Phi in the carrier produces the same phase change in the beat signal.} 

^{-------------------}

^{matlab code from GPT-4 to look at phase variations:} 

%% HFcoherence.m

% C. Stubbs Aug 3 with GPT-4 help

% Read in the WAV file

infile='/Users/cstubbs/data/ionosphere/RMV.wav';

[signal, Fs] = audioread(infile);

% Parameters

Fs=12000; % sample rate in Hz

beat=500; % guess of beat freq

Fn = Fs/2;                                                  % Nyquist Frequency

T = 1/Fs; % Sampling period in seconds

L = length(signal); % Length of the signal, in samples

t = (0:L-1)*T; % Time vector

Ts = mean(diff(t));      

PointsPerCycle=floor(Fs/beat)

windowSize = PointsPerCycle*20; % Size of the sliding window, make it 20 periods

overlap = floor(windowSize/2);    % Overlap between successive windows

stepSize = windowSize - overlap; % Step size for sliding window

Fv = linspace(0, 1, fix(windowSize/2)+1)*Fn; 

% Preallocate arrays to store amplitude and phase

numWindows = floor((length(signal) - overlap) / stepSize);

amplitude = zeros(1, numWindows);

phase = zeros(1, numWindows);

maxfreq=phase;

% Iterate through the signal with a sliding window

for w = 1:numWindows

    % Define the current window

    startIdx = (w - 1) * stepSize + 1;

    endIdx = startIdx + windowSize - 1;

    window = signal(startIdx:endIdx);

    window=window-mean(window);

    % Compute the Fourier transform of the window

    fftWindow = fft(window);

    % Find the dominant component

    [~, idx] = max(abs(fftWindow));

    idx=18;

    maxfreq(w)=Fv(idx);

    % Compute amplitude and phase of the dominant component

    amplitude(w) = abs(fftWindow(idx));

    phase(w) = angle(fftWindow(idx));

    figure(1)

    plot(Fv(2:end),abs(fftWindow(1:end/2)),'ko-')

    xlim([0 1000])

    grid on

    pause(0.01)

    shg

end

%% Plot amplitude and phase

figure(5)

plot(maxfreq,'ko-')

grid on

shg

figure(10);

subplot(2,1,1);

plot(amplitude);

title('Amplitude of Dominant Fourier Component');

xlabel('Window');

ylabel('Amplitude');

grid on

subplot(2,1,2);

plot(phase);

title('Phase of Dominant Fourier Component');

xlabel('Window');

ylabel('Phase (radians)');

grid on

%%

figure(20)

comet(amplitude,phase)

shg

figure(30)

plot(amplitude,phase,'ko')

grid on

shg

----------------------------

This does a sliding FFT down the data stream, which is 3 minutes long. 
480 data points per FFT window, each of duration 0.377 sec. Adjacent data points are correlated, in triplets due to window overlaps.  
There are about 9000 windows over the 3 minute data file. 

I computed phase within each window at a fixed frequency of 425 Hz. If this is wrong it just produces a steady drift in phase, can fit that out. 

Used the unwrap() function in MATLAB to generate evolution of phase, correcting for 2pi rollover. Did a fit to a straight line and computed residuals from that. 
evolution of phase deviations, y axis is cycles. X axis spans duration of 3 min: