%例程6－6  計算雷達回波信號的多譜勒時頻譜
% e.g.6-6.m for example6-6;
N=512; 
 [fm,am,iflaw]=doppler(N,200,65,10,50); 
 figure;
 plot(real(am.*fm));  
 title( 'Signal in time')
 xlabel('Time');
 ylabel('real part');
 figure;
 
 plot(iflaw); 
 title( 'Instantaneous frequency estiwmation')
 xlabel('Time');
 ylabel('Nomized Frequency (Hz)');
 
 sig=sigmerge(am.*fm,noisecg(N),10);
 figure;
 plot(real(sig));
 title( 'Signal in time')
 xlabel('Time');
 ylabel('real part');
 
 [ifhat,t]=instfreq(sig,11:502,10);
 figure;
 hold on; plot(t,ifhat,'g')
 title( 'Instantaneous frequency estiwmation')
 xlabel('Time');
 ylabel('Nomized Frequency (Hz)');
 
 
 function [fm,am,iflaw]=doppler(N,Fs,f0,d,v,t0,c);
%DOPPLER Generate complex Doppler signal.
%    [FM,AM,IFLAW]=DOPPLER(N,FS,F0,D,V,T0,C) 
%    Returns the frequency modulation (FM), the amplitude 
%    modulation (AM) and the instantaneous frequency law (IFLAW) 
%    of the signal received by a fixed observer from a moving target 
%    emitting a pure frequency f0.
%
%    N  : number of points.  
%    FS : sampling frequency (in Hertz).  
%    F0 : target   frequency (in Hertz).  
%    D  : distance from the line to the observer (in meters).  
%    V  : target velocity    (in m/s)
%    T0 : time center                  (default : N/2).  
%    C  : wave velocity      (in m/s)  (default : 340). 
%    FM : Output frequency modulation.  
%    AM : Output amplitude modulation.  
%    IFLAW : Output instantaneous frequency law.
% 
%   Example: 
%    N=512; [fm,am,iflaw]=doppler(N,200,65,10,50); 
%    subplot(211); plot(real(am.*fm)); 
%    subplot(212); plot(iflaw);
%        [ifhat,t]=instfreq(sigmerge(am.*fm,noisecg(N),15),11:502,10);
%        hold on; plot(t,ifhat,'g'); hold off;
%
%   See also DOPNOISE 
 
%   F. Auger, July 94, August 95 - O. Lemoine, October 95.
%   Copyright (c) 1996 by CNRS (France).
%
%   ------------------- CONFIDENTIAL PROGRAM -------------------- 
%   This program can not be used without the authorization of its
%   author(s). For any comment or bug report, please send e-mail to 
%   f.auger@ieee.org 
 
if (nargin <= 4),
 error ( 'At least 5 parameters are required' ); 
elseif (nargin == 5),
 t0=N/2; c=340.0;
elseif (nargin == 6),
 c=340.0;
end;
 
if (N <= 0),
 error ('The signal length N must be strictly positive' );
elseif (d <= 0.0),
 error ('The distance D must be positive' );
elseif (Fs < 0.0),
 error ('The sampling frequency FS must be positive' );
elseif (t0<1) | (t0>N),
 error ('T0 must be between 1 and N');
elseif (f0<0)|(f0>Fs/2),
 error ('F0 must be between 0 and FS/2');
elseif (v<0),
 error ('V must be positive');
else
 tmt0=((1:N)'-t0)/Fs;
 dist=sqrt(d^2+(v*tmt0).^2);
 fm = exp(j*2.0*pi*f0*(tmt0-dist/c));
 if (nargout>=2), 
  if abs(f0)<eps,
   am=0;
  else
   am= 1.0 ./ sqrt(dist); 
  end
 end;
 if (nargout==3), iflaw=(1-v^2*tmt0./dist/c)*f0/Fs; end;
end ;
                   
function noise=noisecg(N,a1,a2)
%NOISECG Analytic complex gaussian noise.
%   NOISE=NOISECG(N,A1,A2) computes an analytic complex gaussian
%   noise of length N with mean 0.0 and variance 1.0. 
%
%   NOISE=NOISECG(N) yields a complex white gaussian noise.
%
%   NOISE=NOISECG(N,A1) yields a complex colored gaussian noise
%   obtained by filtering a white gaussian noise through a
%       sqrt(1-A1^2)/(1-A1*z^(-1)) 
%   first order filter.
%
%   NOISE=NOISECG(N,A1,A2) yields a complex colored gaussian noise
%   obtained by filtering a white gaussian noise through a
%       sqrt(1-A1^2-A2^2)/(1-A1*z^(-1)-A2*z^(-2)) 
%   second order filter.
%
%   Example :
%    N=512;noise=noisecg(N);mean(noise),std(noise).^2
%    subplot(211); plot(real(noise)); axis([1 N -3 3]);
%    subplot(212); f=linspace(-0.5,0.5,N); 
%    plot(f,abs(fftshift(fft(noise))).^2);
%       
%   See also RAND, RANDN, NOISECU.
 
%   O. Lemoine, June 95/May 96 - F. Auger, August 95.
%   Copyright (c) 1996 by CNRS (France).
%
%   ------------------- CONFIDENTIAL PROGRAM -------------------- 
%   This program can not be used without the authorization of its
%   author(s). For any comment or bug report, please send e-mail to 
%   f.auger@ieee.org 
 
if (N <= 0),
 error ('The signal length N must be strictly positive' );
end;
 
if (nargin==1),
 if N<=2,
  noise=(randn(N,1)+j*randn(N,1))/sqrt(2); 
 else
  noise=randn(2^nextpow2(N),1); 
 end
elseif (nargin==2),
 if (abs(a1)>=1.0),
  error ('for a first order filter, abs(a1) must be strictly lower than 1');
 elseif (abs(a1)<=eps),
  if N<=2,
   noise=(randn(N,1)+j*randn(N,1))/sqrt(2); 
  else
   noise=randn(N,1);
  end
 else
  if N<=2,
   noise=(randn(N,1)+j*randn(N,1))/sqrt(2); 
  else
   Nnoise=ceil(N-2.0/log(a1));
   noise=randn(2^nextpow2(Nnoise),1);
  end
  noise=filter(sqrt(1.0-a1^2), [1 -a1],noise);
 end;
elseif (nargin==3),
 if any(roots([1 -a1 -a2])>1),
  error('unstable filter');
 else
  if N<=2,
   noise=(randn(N,1)+j*randn(N,1))/sqrt(2); 
  else
   Nnoise=ceil(N-2.0/log(max(roots([1 -a1 -a2]))));
   noise=randn(2^nextpow2(Nnoise),1);
  end
  noise=filter(sqrt(1.0-a1^2-a2^2), [1 -a1 -a2],noise);
 end;
end;
 
if N>2,
 noise=hilbert(noise)/std(noise)/sqrt(2);
 noise=noise(length(noise)-(N-1:-1:0));
end
 

function sig=sigmerge(x1,x2,ratio);
%SIGMERGE Add two signals with given energy ratio in dB.
%   SIG=SIGMERGE(X1,X2,RATIO) adds two signals so that a given
%   energy ratio expressed in deciBels is satisfied.
% 
%   X1, X2 : input signals.
%   RATIO  : Energy ratio in deciBels   (default : 0 dB).
%   X      : output signal.
%   X= X1+H*X2, such that 10*log(Energy(X1)/Energy(H*X2))=RATIO
%
%   Example : 
%    sig=fmlin(64,0.01,0.05,1); noise=hilbert(randn(64,1));
%    SNR=15; x=sigmerge(sig,noise,SNR);
%    Esig=mean(abs(sig).^2); Enoise=mean(abs(x-sig).^2);
%    10*log10(Esig/Enoise)
 
%   F. Auger, July 1995.
%   Copyright (c) 1996 by CNRS (France).
%
%   ------------------- CONFIDENTIAL PROGRAM -------------------- 
%   This program can not be used without the authorization of its
%   author(s). For any comment or bug report, please send e-mail to 
%   f.auger@ieee.org 
 
if (nargin<2)
 error('At least two parameters are required');
elseif nargin==2,
 ratio=0;
end;
 
[x1row,x1col] = size(x1);
[x2row,x2col] = size(x2);
 
if (x1col~=1)|(x2col~=1),
 error('X1 and X2 must have only one column');
elseif (x1row~=x2row),
 error('X1 and X2 must have the same number of rows');
elseif (length(ratio)~=1),
 error('RATIO must be a scalar');
elseif (ratio==inf),
 sig = x1;
else
 Ex1=mean(abs(x1).^2);
 Ex2=mean(abs(x2).^2);
 h=sqrt(Ex1/(Ex2*10^(ratio/10)));
 sig=x1+h*x2;
end;