examsigw.m

阵列信号处理的工具箱

function [SNR] = examsigw(sigIn, rangeIx, wavefEnds, channel, pulse, procFlags, printFlag, figFileNames, figSize)

%EXAMSIGW Examines a received radar signal waveform in various ways.
%
%--------
%Synopsis:
%  examsigw(sigIn, rangeIx, wavefEnds)
%  examsigw(sigIn, rangeIx, wavefEnds, channel, pulse, procFlags, printFlag,
%    figFileNames)
%  SNR = examsigw(sigIn, rangeIx, wavefEnds, channel, pulse, procFlags,
%    printFlag, figFileNames)
%
%Description:
%  Examines a received radar signal waveform in various ways and plots graphs.
%  The frequency spectrum, amplitude, phase and instantaneous frequency
%  can be examined.
%
% Possible processing and plotting:
%   (1)  Estimate SNR.
%   (2)  Frequency spectrum.
%   (3)  Mean value of amplitude.
%   (4)  Standard deviation of amplitude.
%   (5)  Phases of all pulses.
%   (6)  Mean value of phase.
%   (7)  Standard deviation of phase.
%   (8)  Differences between start phases.
%   (9)  Difference between start phase and end phase of previous pulse.
%   (10) Estimation of instantaneous frequency. A single pulse.
%   (11) Mean value of instantaneous frequency.
%   (12) Standard deviation of instantaneous frequency.
% Which processings to perform are indivually choosen by the input parameter
%   "procFlags". See below.
%
% (1)  Estimate SNR.
%   The ranges 81:256 are assumed to contain only noise.
%   The ranges 11:73 are assumed to contain only target signal.
%   The SNR is estimated as as the ratio of the mean of the power of the
%   target and the mean of the power of the noise in the above range gates.
%
%Output and Input:
%  sigIn (RxRadarSigT): Signal to examine.
%  rangeIx (IndexT): What range bins to plot.
%  wavefEnds (IndexT): A two element vector with range bin indices for start
%    and end of the waveform to examine.
%  channel (IntScalarT):
%  pulse (IntScalarT):
%  procFlags (Vector of BoolT): An integer vector where each element can
%    be eather true (=1) or false (=0). True means that the corresponding
%    processing (se Description above) shall be done.
%  printFlag (BoolT): True (=1) if to print the plotted graphs to files,
%    else false (=0).
%  figFileNames (StringT): The beginning of the file name when printing
%    the plotted graphs to files.
%
%--------
%Notations:
%  Data type names are shown in parentheses and they start with a capital
%  letter and end with a capital T. Data type definitions can be found in [1]
%  or by "help dbtdata".
%  [D] = This parameter can be omitted and then a default value is used.
%  When the [D]-input parameter is not the last used in the call, it must be
%  given the value [], i.e. an empty matrix.
%  ... = There can be more parameters. They are explained under respective
%  metod or choice.
%
%Examples:
%  examsigw(sig, 9:76, [11 74], 6, 20 ,[0 1 zeros(1,10)],1,'specGraph')
%    % examsigw(sigIn, rangeIx, wavefEnds, channel, pulse, procFlags,
%    % printFlag, figFileNames, figSize)
%
%  %See also reference [2] for examples of the use of this function.
%
%Software Quality:
%  (About what is done to ascertain software quality. What tests are done.)
%
%Known Bugs:
%
%References:
%  [1]: Bj鰎klund S.: "DBT, A MATLAB Toolbox for Radar Signal Processing.
%    Reference Guide", FOA-D--9x-00xxx-408--SE, To be published.
%  [2]: Bj鰎klund S.: "Inledande radarliknande m鋞ningar med FOA:s
%    experimentella digitala gruppantenn", FOA-D--9x-00xxx-408--SE,
%    To be published.
%
%See Also:
%  sigplot2, expsig1, simradarsig, compsim4, compsim5

%   *  DBT, A Matlab Toolbox for Radar Signal Processing  *
% (c) FOA 1994-2000. See the file dbtright.m for copyright notice.
%
%  Start        : 980710 Svante Bj鰎klund (svabj).
%  Latest change: $Date: 2000/11/17 12:45:12 $ $Author: svabj $.
%  $Revision: 1.5 $
% *****************************************************************************


% ----------------------------------------------------------------------- %
% Handle input parameters
% ----------------------------------------------------------------------- %

if (nargin < 3)
  error('DBT-Error: To few input parameters.')
end

% ****************** Add missing input parameters ******************

arginNo=4;
if (nargin < arginNo)
  channel = [];
end
arginNo = arginNo +1;
if (nargin < arginNo)
  pulse = [];
end
arginNo = arginNo +1;
if (nargin < arginNo)
  procFlags = [];
end
arginNo = arginNo +1;
if (nargin < arginNo)
  printFlag = [];
end
arginNo = arginNo +1;
if (nargin < arginNo)
  figFileNames = [];
end
arginNo = arginNo +1;
if (nargin < arginNo)
  figSize = [];
end
arginNo = arginNo +1;


% ****************** Default values ******************

if isempty(channel)
  channel = 1;
end%if
if isempty(pulse)
  pulse = 1;
end%if
if isempty(procFlags)
  procFlags = ones(1,12);
end%if
if isempty(printFlag)
  printFlag = 0;
end%if
if isempty(figFileNames)
  figFileNames = 'examsigw';
end%if
if isempty(figSize)
  figSize = 8;
end%if

% ****************** Error check input parameters ******************


% ----------------------------------------------------------------------- %
% Parameters.
% ----------------------------------------------------------------------- %

%sig11 = sig;
% With this (uncompensated) signal there is small difference of the amplitude levels but the amplitude shapes look the same.

%sig11 = sigIn;
% With this (compensated) signal there is not notable difference of the amplitude levels or shapes.


%rangeIx = 9:76;
%wavefEnds = [3 66];
waveIx = wavefEnds(1):wavefEnds(2);


%printFlag = 0		% 1 = print to eps-files. 0 = do not print.

%channel = 6;
set(0,'DefaultAxesFontSize',8)
set(0,'DefaultAxesFontAngle','italic')

plotStartShift = -wavefEnds(1);
plotIx = rangeIx + plotStartShift;
plotDiffIx = plotIx(1:length(plotIx)-1)+1;
plotWaveEnds = wavefEnds + plotStartShift;
plotDiffWaveEnds = [plotWaveEnds(1), plotWaveEnds(2)-1];

%plotStartIx = -1;
%plotStartShift = plotIx + plotStartIx;


% ----------------------------------------------------------------------- %
% Estimate SNR.
% ----------------------------------------------------------------------- %

procNo=1;
if procFlags(procNo)
  %pulse = 20;
  %channel = 6;
  noisePower = mean(abs(sigIn.signals(pulse,81:256,channel)).^2);
    % The ranges 81:256 are assumed to contain only noise.
  tgtPower = mean(abs(sigIn.signals(pulse,11:73,channel)).^2);
    % The ranges 11:73 are assumed to contain only target signal.
  SNR = 10*log10(tgtPower/noisePower)
end%if procFlags



% ----------------------------------------------------------------------- %
% Plot frequency spectrum.
% ----------------------------------------------------------------------- %

procNo=2;
if procFlags(procNo)
  specLen = 2048;
  fSamp =  1/sigIn.waveform.sampleTime;
  %spect = fftshift(abs(fft((sigVec))));
  %spect = fftshift(abs(fft(sigIn.signals(10,1:256,channel))));
  spect = fftshift(abs(fft(sigIn.signals(10,waveIx,channel),specLen)));

  absSpec = abs(spect);
  f = linspace(-fSamp/2,fSamp/2,specLen);
  figure,plot(f,absSpec/max(absSpec))
  plotxline([-fSamp/2, 0, fSamp/2],'r','--')
  title('Frequency spectrum')
  xlabel('Frequency [Hz]')
  ylabel('Normalized spectrum')
  if (printFlag)
    printfm( [figFileNames, num2str(procNo)],[], figSize)
  end%if
end%if procFlags

% ----------------------------------------------------------------------- %
% Statistical analysis.
% ----------------------------------------------------------------------- %

sSize = sigsize(sigIn);
noPulses = sSize(1);
noRanges = length(rangeIx);
freqEstMat = zeros(noPulses, noRanges-1);
sigAmpMat = zeros(noPulses, noRanges);
sigPhaseMat = zeros(noPulses, noRanges);
for pulse = 1:noPulses
    sigVec = sigIn.signals(pulse,rangeIx,channel);
    sigAmpMat(pulse,:) = abs(sigVec);
    sigPhaseMat(pulse,:) = unwrap(angle(sigVec));
    freqEst = estlocalfreq(sigVec, fSamp);
    freqEst = freqEst(:).';
    freqEstMat(pulse,:) = freqEst;
end%for pulse
sigAmpMean = mean(sigAmpMat);
sigAmpStd = std(sigAmpMat);
sigPhaseMean = mean(sigPhaseMat);
%startPhases = sigPhaseMat(:,wavefEnds(1));
startPhases = sigPhaseMat(:,3);			% NOTE: Should use a variable
						% instead of number constant.
startPhasesUW = unwrap(startPhases);
startPhasesDiff = diff(startPhasesUW);
%r2d(startPhases).'
%r2d(startPhasesDiff).'

%'a'
%wavefEnds(1)
%wavefEnds(2)
%size(sigPhaseMat)
endPhases = sigPhaseMat(:,66);			% NOTE: Should use a variable
						% instead of number constant.
endPhasesUW = unwrap(endPhases);
%diffStartEndPhases = startPhasesUW(2:noPulses) - endPhasesUW(1:noPulses-1);
startRangeIx = 3;
endRangeIx = 66;
diffStartEndPhases = unwrap(angle(sigIn.signals(2:noPulses,startRangeIx,channel) ./ sigIn.signals(1:noPulses-1,endRangeIx,channel)));
%r2d(diffStartEndPhases)
%endPhasesDiff = diff(endPhasesUW);

sigPhaseMat = sigPhaseMat - repmat(startPhases,1 , size(sigPhaseMat,2));
sigPhaseStd = std(sigPhaseMat);
freqMean = mean(freqEstMat);
freqStd = std(freqEstMat);

maxSigAmpMean = max(sigAmpMean);
%size(maxSigAmpMean)

% ----------------------------------------------------------------------- %
% Plot mean value of amplitude.
% ----------------------------------------------------------------------- %

procNo=3;
if procFlags(procNo)
  figure
  plot(plotIx,sigAmpMean/maxSigAmpMean)
  xlabel('Time in pulse')
  title('Amplitude (mean)')
  ylabel('Normalized amplitude')
  %ylim([0 0.6e6])
  plotxline(plotWaveEnds,'r','--')
  % Typical look for several channels.
  %if (printFlag)
  %  printfm m1079aproc3.eps [] 8
  % end%if
  if (printFlag)
    printfm( [figFileNames, num2str(procNo)],[], figSize)
  end%if
end%if procFlags

% ----------------------------------------------------------------------- %
% Plot standard deviation of amplitude.
% ----------------------------------------------------------------------- %

procNo=4;
if procFlags(procNo)
  figure
  plot(plotIx,sigAmpStd/maxSigAmpMean)
  xlabel('Time in pulse')
  title('Amplitude (standard deviation)')
  ylabel('Normalized amplitude')
  %ylim([0 0.6e6])
  plotxline(plotWaveEnds,'r','--')
  %if (printFlag)
  %  printfm m1079aproc4.eps [] 8
  %end%if
  if (printFlag)
    printfm( [figFileNames, num2str(procNo)],[], figSize)
  end%if
end%if procFlags



% ----------------------------------------------------------------------- %
% Plot phases of all pulses.
% ----------------------------------------------------------------------- %

procNo=5;
if procFlags(procNo)
  figure
  plot(plotIx,r2d(sigPhaseMat.'))
  xlabel('Time in pulse')
  title('Phases of all pulses. Normalized start phase.')
  ylabel('[degrees]')
  %ylim([0 0.6e6])
  plotxline(plotWaveEnds,'r','--')
  if (printFlag)
    printfm( [figFileNames, num2str(procNo)],[], figSize)
  end%if
end%if procFlags



% ----------------------------------------------------------------------- %
% Plot mean value of phase.
% ----------------------------------------------------------------------- %

procNo=6;
if procFlags(procNo)
  figure
  plot(plotIx,r2d(sigPhaseMean))
  xlabel('Time in pulse')
  title('Phase (mean)')
  ylabel('[degrees]')
  %ylim([0 0.6e6])
  plotxline(plotWaveEnds,'r','--')
  if (printFlag)
    printfm( [figFileNames, num2str(procNo)],[], figSize)
  end%if
end%if procFlags


% ----------------------------------------------------------------------- %
% Plot standard deviation of phase.
% ----------------------------------------------------------------------- %

procNo=7;
if procFlags(procNo)
  figure
  plot(plotIx,r2d(sigPhaseStd))
  xlabel('Time in pulse')
  title('Phase (standard deviation)')
  ylabel('[degrees]')
  ymax(2)
  %ylim([0 0.6e6])
  plotxline(plotWaveEnds,'r','--')
  if (printFlag)
    printfm( [figFileNames, num2str(procNo)],[], figSize)
  end%if
end%if procFlags


% ----------------------------------------------------------------------- %
% Plot differences between start phases.
% ----------------------------------------------------------------------- %

procNo=8;
if procFlags(procNo)
  figure
  plot(1:127,r2d(startPhasesDiff.'))
  xlabel('Pulses')
  title('Diff Start Phases')
  ylabel('[degrees]')
  %ylim([0 0.6e6])
  plotxline(plotDiffWaveEnds,'r','--')

  fprintf('Mean of start phase diff %g degrees.\n', ...
  r2d(mean(startPhasesDiff)))
  fprintf('Standard deviation of start phase diff %g degrees.\n', ...
  r2d(std(startPhasesDiff)))
  if (printFlag)
    printfm( [figFileNames, num2str(procNo)],[], figSize)
  end%if
end%if procFlags


% ----------------------------------------------------------------------- %
% % Plot difference between start phase and end phase of previous pulse.
% ----------------------------------------------------------------------- %

procNo=9;
if procFlags(procNo)
  figure
  plot(1:127,r2d(diffStartEndPhases.'))
  xlabel('Pulses')
  title('Diff Start End Phases')
  ylabel('[degrees]')
  %ylim([0 0.6e6])
  plotxline(plotDiffWaveEnds,'r','--')
  if (printFlag)
    printfm( [figFileNames, num2str(procNo)],[], figSize)
  end%if
end%if procFlags


% ----------------------------------------------------------------------- %
% Plot estimation of instantaneous frequency. A single pulse.
% ----------------------------------------------------------------------- %

procNo=10;
if procFlags(procNo)
  sigVec = sigIn.signals(10,rangeIx,channel);
  freqEst = estlocalfreq(sigVec, fSamp);
  freqEst = freqEst(:).';
  len = length(sigVec);
  Ts = sigIn.waveform.sampleTime;
  t = linspace(0, Ts*len, len);
  figure
  %plot(t,[0,freqEst])
  %xlabel('Time [s]')
  plot(plotDiffIx,freqEst/1000)
  xlabel('Time in pulse')
  title('Frequency estimation')
  ylabel('Frequency [kHz]')
  ylim([-50 600])
  %ylim([0 0.6e6])
  plotxline(plotDiffWaveEnds,'r','--')
  if (printFlag)
    printfm( [figFileNames, num2str(procNo)],[], figSize)
  end%if
end%if procFlags


% ----------------------------------------------------------------------- %
% Plot mean value of instantaneous frequency.
% ----------------------------------------------------------------------- %

procNo=11;
if procFlags(procNo)
  figure
  plot(plotDiffIx,freqMean/1000)
  xlabel('Time in pulse')
  title('Frequency estimation (mean)')
  ylabel('Frequency [kHz]')
  %ylim([0 0.6e6])
  ylim([-50 600])
  plotxline(plotDiffWaveEnds,'r','--')
  % Typical look for several channels.
  %if (printFlag)
  %  printfm m1079aproc9.eps [] 8
  %end%if
  if (printFlag)
    printfm( [figFileNames, num2str(procNo)],[], figSize)
  end%if
end%if procFlags

% ----------------------------------------------------------------------- %
% Plot standard deviation of instantaneous frequency.
% ----------------------------------------------------------------------- %

procNo=12;
if procFlags(procNo)
  figure
  plot(plotDiffIx,freqStd/1000)
  xlabel('Time in pulse')
  title('Frequency estimation (standard deviation)')
  ylabel('Frequency [kHz]')
  %ylim([0 1e4])
  ylim([0 7])
  plotxline(plotDiffWaveEnds,'r','--')
  % Typical look for several channels.
  if (printFlag)
    printfm( [figFileNames, num2str(procNo)],[], figSize)
  end%if
end%if procFlags