format short
format compact
echo on
clc
% Let's create a signal consisting of 1024 samples
% of a frequency-modulated sinusoid,
% ramping linearly from 128 Hz, (assuming 1024 samples per second)
% to 256 Hz.
%
Nsamps = 1024;
Fstart = 128;
Fdiff = 128;
F = Fstart + Fdiff*(0:(Nsamps-1))/Nsamps;
A = ones(1,Nsamps);
S1 = afmod(A,F);
figure(1)
subplot(4,1,1)
plot(A); title('Amplitude')
subplot(4,1,2)
plot(F); title('Frequency')
subplot(4,1,3)
plot(S1(1:200)); title('First 200 samples of F modulated sinusoid')
subplot(4,1,4)
plot(S1(824:1024)); title('Last 200 samples of F modulated sinusoid')
pause
clc
% If we use an DFT of the whole sequence
% to show us the frequency content of this signal,
% its essential characteristics will not be very plainly shown.
% An amplitude spectrum 
%    tells us what range of frequencies were involved:
Spec = fft(S1);
figure(2);
subplot(2,1,1)
plot(abs(Spec(1:(Nsamps/2))));
title('Amplitude Spectrum: linearly F-modulated Sinusoid');
% The phase spectrum provides interesting information,
% but it certainly doesn't give us a direct picture:
subplot(2,1,2)
plot(unwrap(angle(Spec(1:(Nsamps/2)))));
title('Phase Spectrum: linearly F-modulated Sinusoid');
pause;
clc
% One way to see that we are missing something 
%   is to compare the amplitude and phase spectra 
%   of a time-reversed version of the signal:
Spec1 = fft(S1(1024:-1:1));
figure(3);
subplot(2,1,1)
plot(abs(Spec1(1:(Nsamps/2))));
title('Amplitude Spectrum: time-reversed signal');
subplot(2,1,2)
plot(unwrap(angle(Spec1(1:(Nsamps/2)))));
title('Phase Spectrum: time-reversed signal');
pause;
clc
% Now let's make a sinusoid whose frequency is modulated sinusoidally.
Fcarrier = 256;
Fmod = 128;
F2 = Fcarrier + makesin(Nsamps,2,Fmod,0);
S2 = afmod(A,F2);
figure(1)
subplot(3,1,1)
plot(F2); title('Sinusoidal frequency modulation: frequency')
subplot(3,1,2)
plot(S2(1:200)); title('First 200 samples of modulated signal')
subplot(3,1,3)
plot(S2(201:400)); title('Second 200 samples')
pause
clc
% Again, the overall DFT is not especially revealing:
Spec2 = fft(S2);
figure(2);
subplot(2,1,1)
plot(abs(Spec2(1:(Nsamps/2))));
title('Amplitude Spectrum of Sinusoidal F-modulated signal');
subplot(2,1,2)
plot(unwrap(angle(Spec2(1:(Nsamps/2)))));
title('Phase Spectrum');
pause;
clc
% Often, a spectrogram-style analysis, in which we make a series of
% time-windowed Dfts, will show pretty well what is going on.
% But when modulation is rapid enough, a time-windowed
% DFT that is wide enough in time to get reasonable frequency resolution
% will contain enough change in frequency to blur the peak considerably:
figure(1)
subplot(2,1,1)
Slocal = S2(128:383);
LocalSpec = fft(hamming(256)'.*Slocal);
plot(512*(1:128)/128,abs(LocalSpec(1:128)))
title('Amplitude spectrum - 256 samples of F-modulated sinusoid')
subplot(2,1,2)
sin256 = makesin(1024,256,1,0);
LS256 = fft(hamming(256)'.*sin256(128:383));
plot(512*(1:128)/128, abs(LS256(1:128)))
title('Amplitude spectrum - 256 samples of unmodulated sinusoid')
hold off;
pause;
clc
%% Let's stop for a minute and think about how we might
%  analyze such an FM/AM sinusoid more directly,
%  i.e. in such a way as to recover the FM and AM time functions directly.
%  In other words, we need a way to estimate the instantaneous
%  frequency and amplitude of a sinusoid.
%
%  Suppose for a moment that we have a sine and cosine of the same
%  frequency,
%  combined as cos + i*sin
%  Then sequential samples of the combined signal
%  will move around the unit circle at a rate 
%  that is equal to their (shared) frequency.
pause
% We can get the frequency back that we put in
%    from two adjacent samples of each signal --
% by calculating the rate of change of the angle
SINE = makesin(1000, 3.7, 1.3, 0);
COS = makecos(1000, 3.7, 1.3, 0);
F = angle(COS(2) + i*SINE(2)) - angle(COS(1) + i*SINE(1))
% and then multiplying by Nsamples/(2*pi) 
%     to get the units right:
F*1000/(2*pi)
pause
% We can get the amplitude back from any pair of samples of the signal:
A = abs(COS(2) + i*SINE(2))
clc
%
% Note that it was essential to look at the rate of change of the
% angle -- if we only knew the sine and cosine values at one point,
% we know the phase, but the frequency could have been anything.
%
% Note for future reference that if the amplitude of one of the
% sinusoids is modified a bit, it will affect not only the estimate
% of amplitude but also the estimate of frequency:

SINE = makesin(1000, 3.7, 1.3, 0);
COS = makecos(1000, 3.7, 1.301, 0);
F = diff(unwrap(angle(COS + i*SINE)))*1000/(2*pi);
F(1:10)'
pause
A = abs(COS + i*SINE);
A(1:10)'
figure(1)
plot(F)
title('Frequency estimates for amp(COS) = 1.301, amp(SINE) = 1.3')
figure(2)
plot(A)
title('Amplitude estimates for amp(COS) = 1.301, amp(SINE) = 1.3')
%
% What if we are given an input sinusoid of unknown frequency,
% amplitude and phase? How might we estimate these quantities?
% Since the sinusoid is assumed to have constant parameters
% we could recover the frequency, amplitude and phase from a DFT --
% we can get increasing accuracy from increasing the number of samples.
% But let's try to find a way of doing it that would generalize to
% to the FM/AM case.
pause
clc
%
% The cosine function is just the sine function phase-shifted by 90 degrees.
% Adding a constant phase shift to both functions would not change
% the derivative of the angle, which is all that we care about in
% reconstructing the frequency, nor would it change the magnitude,
% which is all the matters in recovering the amplitude.
% So to use the same technique as before, all we need is to get
% a version of our input sinusoid that is shifted in phase by 90 degrees
% at all frequencies.
% Since phase shifting is a linear operation, there is some linear
% operator that does exactly this.
% It is an important enough operator to have a special name:
%
%   Hilbert transformer
%
pause
clc
% 
% What is this operator?
%
% There are several ways to approach it.
% For instance, we could take a DFT, 
%     adjust the phase in the frequency domain,
%     and then take the inverse DFT.
% This is what the Matlab function hilbert() does -- here's how:
%
% "It turns out that"
% if a sequence is CAUSAL (i.e. its left half is zero)
% then the real and imaginary parts of its Fourier transform
%      are in the Hilbert transform relation.
%
% Thus we can create a Hilbert transform pair by
%
%   (1) taking an FFT
%   (2) zeroing out the negative frequencies
%   (3) taking an inverse FFT
%
% This should produce a complex sequence whose real part is the original,
% and whose imaginary part is the original phase-shifted by 90 degrees.
%
%
figure(1)
subplot(1,1,1)
hCOS = hilbert(COS);
plot(real(hCOS)); title('Real & imaginary parts of hilbert(COS)')
hold on
plot(imag(hCOS),':')
hold off
%
% As we can see, it approximately works, 
%    although there are some significant edge effects.
pause
% In the middle, the estimates are not too bad, 
%    but on the edges, they are not so good:

FhCOS = diff(unwrap(angle(hCOS)))*1000/(2*pi);
FhCOS(501)
FhCOS(2)
abs(hCOS(501))
figure(2)
subplot(2,1,1)
plot(FhCOS); title('Frequency 3.7 recovered by hilbert()')
subplot(2,1,2)
plot(abs(hCOS)); title('Amp 1.301 recovered by hilbert()')
pause
clc
% Aside from the edge effects, which we might try to
% mitigate by windowing or by other means,
% this approach has got the difficulty that we have to DFT the entire sequence.
% Since we might be interested in arbitrarily-long sequences,
% this could get tedious.
% Since the Hibert transformer is a linear shift-invariant operation,
% there must be a convolutional operator that accomplishes it.
%
pause
% This operator is called a "Hilbert transformer".
% "It turns out that" 
%     the impulse response of a Hilbert transformer is:
%
% (2/pi)* ( sin(pi*n/2).^2 )/n
%
% where -inf < n < inf
%
% Because this extends to infinity in both directions,
% the ideal Hilbert transformer, like the ideal low-pass filter,
% is not realizable.
%
% We have to make some compromise, either by windowing the impulse response
% or by some other method.
pause
clc
% A windowed Hilbert impulse response will work pretty well,
% as long as the window is larger than the local period of the sinusoid:
L = 401;
n = (L-1)/2;
m = -n:n;
hh = sin(pi*m/2).^2;
hh = (2/pi)*hh./m;
% sorry about the divide by zero!
hh(n+1) = 0;
h1 = hh.*hanning(L)';
figure(1)
lplot(h1); title('Impulse response of windowed Hilbert Transformer')
pause
clc
% Let's try it out; first on the simple cosine wave:
hCOS = conv(COS,h1);
N = length(h1);
Nsamps = length(COS);
pause
% We need to throw away the extra samples created by convolution,
hCOS = hCOS( ((N+1)/2):(Nsamps+(N-1)/2) );
Fnew = diff(unwrap(angle(i*hCOS + COS))) * Nsamps/(2*pi);
Fnew(501)
figure(2)
subplot(2,1,1)
plot(Fnew); title('Frequency of COS recovered by Hilbert transformer')
pause
% This method will not work so well if the hilbert transformer
% is small relative to the local period:
h1 = s90(17)
hCOS = conv(COS,h1);
N = length(h1);
Nsamps = length(COS);
hCOS = hCOS( ((N+1)/2):(Nsamps+(N-1)/2) );
Fnew = [0 diff(unwrap(angle(i*hCOS + COS))) * Nsamps/(2*pi)];
Fnew(501)
subplot(2,1,2)
plot(Fnew); title('Frequency of COS (not) recovered by too-short Hilbert operator')
hold on
plot(COS,':');
hold off
pause
clc
% Now let's try something more interesting --
% Remember that S1 is the signal with a ramp in frequency...
h1 = s90(17);
N = length(h1);
Sh = conv(S1,h1);
% again we get rid of the extra samples...
% and we might as well multiply by i at the same time...
Nsamps = length(S1);
Shi = i * Sh(((N+1)/2):(Nsamps+(N-1)/2));
Fh = diff(unwrap(angle(Shi+S1)));
figure(1);
subplot(1,1,1)
plot(Fh * Nsamps/(2*pi)); title('Frequency ramp of S1 recovered by Hilbert')
axis([-50 1100 100 280]);
pause
clc
%
% Now let's try S2, where the f-modulation was sinusoidal:
Sh2 = conv(S2,h1);
Shi2 = i * Sh2(((N+1)/2):(Nsamps+(N-1)/2));
Fh2 = diff(unwrap(angle(Shi2 + S2)));
figure(2)
subplot(2,1,1)
plot(Fh2 * Nsamps/(2*pi));title('Frequency of S2 recovered by Hilbert')
subplot(2,1,2)
plot(abs(Shi2 + S2));title('Amplitude of S2 recovered by Hilbert')
pause
clc
% Now let's try it 
%      with some real amplitude modulation added in...
A2 = [(1:100)/100   ones(1,100) ...
      (1 + makesin(624,1,.5,0)) ...
       ones(1,100)   (100:-1:1)/100 ];
S3 = S2 .* A2;
figure(1)
subplot(2,1,1)
plot(F2); title('Frequency modulation of S2')
subplot(2,1,2)
plot(A2); title('Amplitude modulation of S2 to produce S3')
pause
clc
% Now try the analysis:
Sh3 = conv(S3,h1);
Shi3 = i * Sh3(((N+1)/2):(Nsamps+(N-1)/2));
Fh3 = diff(unwrap(angle(Shi3 + S3)));
Ah3 = abs(Shi3+S3);
figure(2)
subplot(2,1,1)
plot(Ah3); title('estimated amplitude modulation');
subplot(2,1,2)
plot(Fh3); title('estimated frequency modulation');
pause
clc
% Now let's try the same analysis
%     in a completely different way,
%     via a recent idea of Jim Kaiser's....
%
% The idea is based on something that he calls
%     Teager's Energy Operator
% which is defined as
%
%    psi[x(n)] = x(n)*x(n) - x(n+1)*x(n-1)
%
pause
% If x is a (sampled) frequency-and-amplitude-modulated sinusoid,
%    with F and A being its the instantaneous frequency and amplitude, then if
%
%         p  = psi(x)
%   and   pd = psi(deriv(x))
%
%   [ where psi(deriv(x)) is (psi(diff(x)) + psi(shift1(diff(x))))/2 ]
%
% then if   a = 1 - (pd ./ (2*p)
%
%    F = acos(a)         A = sqrt(abs(p ./ (1 - a.^2)))
%
pause
clc
[xam xfm] = dsa1(S3);

figure(3)
subplot(2,1,1)
plot(xam(3:1022)); title('Amplitude by Kaiser/Teager')
subplot(2,1,2)
plot(xfm(3:1022));; title('Frequency by Kaiser/Teager')
pause
clc
% Kaiser's method works pretty well regardless of the frequency
% of the sinusoid, as long as there is only one component.