%   Propagate a pulse to simulate the non-linear distortion
%  and show it for pulse inversion
%
%  Calling: pulse = propagate_pulse_inversion (P, z_max, fs)
%
%  Input:  P     - Pulse to propagate
%          z_max - Distance to propagate
%          fs    -  Sampling frequency
%
%   Version 1.1, April 22, 2003, JAJ
%   Version 1.2, October 1, 2018, JAJ
%      Adapted for newer versions of Matlab

function  pulse = display_pulse_inversion (hp, hdft, P, z_max, fs)

%  Propagate the pulse

Pn=propagate_pulse (P, z_max, fs);
Pnn=propagate_pulse (-P, z_max, fs);

%  Display the result

axes(hp)
plot((1:length(Pn))/fs,Pn,'r',(1:length(Pnn))/fs,Pnn,'g',(1:length(Pnn))/fs,Pn+Pnn,'b','LineWidth',2)
xlabel('Relative time [s]')
ylabel('Pressure [Pa]')
title(['Non-linear pulses at ',num2str(z_max*100),' cm'])
legend('Positive non-linear pulse','Negative non-linear pulse','Pulse inversion','Location','SouthEast');

axes(hdft)
dft_log=20*log10(abs(fft(Pn,1024)));
dft_logn=20*log10(abs(fft(Pn+Pnn,1024)));
dft_logn=dft_logn-max(dft_log);
dft_log=dft_log-max(dft_log);
plot((0:1023)/1024*fs/1e6,dft_log,'r--',(0:1023)/1024*fs/1e6,dft_logn,'b','LineWidth',2)
xlabel('Frequency [MHz]')
ylabel('Spectral amplitude [dB]')
axis([0 30 -60 0])
grid
title('Spectra of pulses')
legend('Non-linear spectrum','Pulse inversion spectrum');
