% A simple Matlab demonstration of detection of harmonic signal in additive
% noise. "An input signal" is additive mixture of randomly shifted segment
% of harmonic signal and Gaussian white noise.
% The detection is based on:
% a) "correlation" principle
% b) by usage of convolution and matched filter.
%
% Main goal of this example is demonstartion of usage of specialized m-functions
% (filter function in this example) and demonstartion of relation between
% correlation and convolution.
% Basic work with external files is also demonstarted.
%
% PEEC, M.D. 2007-10-08, v.1.0

fvz = 10000;       % sampling frequency
fsig = 200;        % fequency of "useful" harmonic signal
Asig = 1.0;        % amplitude of "uzitocneho" harmonic signal
Asum = 2.0;        % a parameter determines noise "amplitude" (!!! it is not amplitude)

N = 1024;          % length of the record (number of samples) to be analyzed
T = 100;           % length (number of samples) of the segment of "useful" signal

% generation and visualization of "useful" signal
sig = Asig*sin(2*pi*fsig*(0:T-1)/fvz);  % generation

figure(1);         % display vill be generated in separate window no.1
plot(sig);         % display command
title('Useful signal');               % picture description
xlabel('Discrete time');
ylabel('Amplitude');
grid;              % switch on the grid in dispalyed picture

pause;             % stop processing, waiting for Enter

% define of "random" position of "useful" signal within (N sample long) record
pol = round(1+ rand(1)*(N-T));  % start position
in = [ zeros(1,pol) sig ];      % "place" useful signal into the final record (from the pol+1 position)
in(N) = 0;                      % definitition (force to zero) of all resulting record samples

figure(2);         % display vill be generated in separate window no.2
plot(in);          % display command
title('Randomly shifted useful signal');  % picture description
xlabel('Discrete time');
ylabel('Amplitude');
grid;              % switch on the grid in dispalyed picture

pause;

% generation of noisy signal (additive mixture of signal in + Gaussian white noise)
sum = Asum*randn(1,N); % generation of Gaussian white noise
in_sum = in + sum; % additive mixture

figure(3);         % display vill be generated in separate window no.3
plot(in_sum);      % display command
title('Randomly shifted noisy signal');  % picture description
xlabel('Discrete time');
ylabel('Amplitude');
grid;              % switch on the grid in dispalyed picture

pause;

% Computation of CORRELATION between (vectors) sig and in_sum. The computation is
% performed very inefficiently by using only basic Matlab commands.
% More efficient way of computation will be decribed later during lessons.

out = zeros(1,N-T+1);  % memory allocation, not necessary but can speed up computation
for k=1:(N-T+1)    % for cycle in Matlab
   out(k) = sig*in_sum(k:k+T-1)';   % computation of correlation for shift "k"
end

% Computation by using convolution (by using matched filter)
N_sig = length(sig);
for k=1:N_sig                 % impulse respones of matched filter is mirrored
   h_imp(k) = sig(N_sig-k+1); % version of the signal to be detetcted
end

% filtration (CONVOLUTION!!!) of noisy signal by using matched filter
% function filter is a standard (highly optimized) Matlab one
in_sum2 = [in_sum zeros(1,T)];   % extension of the signal by zero vector appending
out2 = filter(h_imp, 1, in_sum);
out2 = out2(N_sig:length(out2)); % remove a transient response at the beginning


figure(4);         % display vill be generated in separate window no.4
plot(out);         % display command
title('Output of correlation receiver');  % picture description
xlabel('Shift');
ylabel('Amplitude');
grid;              % switch on the grid in dispalyed picture

pause;             % stop processing, waiting for Enter

% picture display with predefined scaling on x nad y axeses
figure(5);         % display vill be generated in separate window no.5
plot(out-out2);    % display command
title('Difference between "convolution" and "correlation" based compuatation');  % picture description
xlabel('Shift');
ylabel('Error');
grid;              % switch on the grid in dispalyed picture

% picture display with pre-defined resolution in x coordinate
% (demonstration of Matlab comand "axis")
roz = axis;        % resolution in x and axes, see Figure(5)

figure(6);         % display vill be generated in separate window no.6
plot(out-out2);    % display command
title('Difference display with defined range in x coordinate');  % picture description
xlabel('Shift');
ylabel('Error');
axis([1 N-T+1 roz(3) roz(4)]);  % roz(3), roz(4) is the resolution in y coordinate
grid;              % switch on the grid in dispalyed picture

pause

[m,odhad] = max(out);  % find a maximum
pol = pol+1;       % real start is shifted by 1 sample
pol                % precise starting position (defined during signals generation)
odhad              % estimation of starting position by processing noisy signal

% demonstration of work with external files. The basic principle is similar as in C-language
% (stdio library). This apporach provides higher flexibility than a very simple Matlab command
% save -ascii.

% Next commands write outut values (after rounding) out as integer numbers
% in standard ASCII format, as well as in a binary format (2 bytes/sample).
% The binary format can used e.g. in direct communication with hardware tools
% (DA converters, microcontrollers, etc.).
x = round( out );   % rounding to integer numbers
save x.txt x -ascii % save samples to the external file in ASCI format (reference computation)
outFile = fopen('x.dat','wb'); % create a file handle (x.dat file)
                               % write data in binary format (wb)

fwrite(outFile,x,'short');     % save samples in binary format (2 bytes/sample)
fclose(outFile);

input('and now compare x.dat and x.txt...');


