Follow the script:
clear all;
close all;
% set the size of the perfectly matched layer (PML)
PML_X_SIZE = 10; % [grid points]
PML_Y_SIZE = 10; % [grid points]
Nx = 1024 ; %
dx = 1e-3; %
Ld = 300; %
kgrid = kWaveGrid(Nx, dx);
posTransm = (Nx-Ld)/2;
posRecep = posTransm + Ld;
medium.sound_speed = 1500 * ones(Nx, 1);
medium.sound_speed(posTransm:posRecep)=4000;
medium.density = 1000 * ones(Nx,1);
medium.density(posTransm:posRecep)= 2200;
t_end = kgrid.x_size / max(medium.sound_speed(:));
dt=75e-9;
Fs=1/dt;
Fmin = 1e5;
Fmax = 5e5;
Lch = 255;
Tfin = Lch*dt;
EixoT = (0:dt:Tfin);
ch = chirp(EixoT, Fmin, Tfin, Fmax);
%Apply Gaussian filter that simule effect of transducer
%ch = gaussianFilter(ch, Fs, 3e5, 10);
%Define source location
source.p_mask=zeros(Nx,1);
source.p_mask(posTransm)=1; %
source.p = ch;
%Define receiver location
sensor.mask = zeros(Nx,1);
sensor.mask(posTransm) = 1;
sensor.mask(posRecep) = 1
%sensor.mask = [-150e-3, 150e-3]; % [mm]
kgrid.t_array = makeTime(kgrid,medium.sound_speed(:),[],t_end);
% Run model
sensor_data = kspaceFirstOrder1D(kgrid, medium, source, sensor, 'PlotLayout', true, 'Smooth', [false, false, false]);
R1=(sensor_data(1, :)); % Output of R1
R2=(sensor_data(2, :)); % Output of R2
% Apply low pass filter
%R1 = applyFilter(R1, 1/dt, 7.5e5, 'LowPass', 'Plot', false, 'ZeroPhase', true);
%R2 = applyFilter(R2, 1/dt, 7.5e5, 'LowPass', 'Plot', false, 'ZeroPhase', true);
disp('Fill signal length with extra zeros up to 2 power closer')
tam=2^(ceil(log2(length(R1)))); %
EixoS = (dt:dt:dt*(tam));%
%R1p = [zeros(length(ch),1);R1';zeros(tam-length(ch)-length(R1),1)];
%R2p = [zeros(length(ch),1);R2';zeros(tam-length(ch)-length(R2),1)];
%%Without moving the beginning of strokes with zeros
R1p = [R1';zeros(tam-length(R1),1)];
R2p = [R2';zeros(tam-length(R2),1)];
chlong = [ch';zeros(tam-(length(ch)),1)];
%chlong = [zeros(length(ch),1);ch1';zeros(tam-2*(length(ch)),1)];
%TF_R1=fft(R1p); % FFT in R1
%TF_R2=fft(R2p); % FFT in R2
%TF_ch = fft(chlong); % FFT in chirp
%R1c = real(ifft(TF_R1/TF_ch)); % Deconvolution R1
%R2c = real(ifft(TF_R2/TF_ch)); % Deconvolution R2
R1c = real(ifft(fft(R1p)./fft(chlong)));
R2c = real(ifft(fft(R2p)./fft(chlong)));
[q1,r1] = deconv(R1,chlong)
% Apply low pass filter
%R1c = applyFilter(R1c, 1/dt, 7.5e5, 'LowPass', 'Plot', false, 'ZeroPhase', true);
%R2c = applyFilter(R2c, 1/dt, 7.5e5, 'LowPass', 'Plot', false, 'ZeroPhase', true);
% -------------------------------------------------------------------------
% Plot results
% -------------------------------------------------------------------------
figure;
subplot(2, 1, 1);
plot (EixoS,R1c,EixoS,R1p);
legend('Signal R1 and Deconvolved');
ylabel('Power');
xlabel('Time (us)')
grid minor;
hold on;
subplot(2, 1, 2);
plot (EixoS,R2c,EixoS,R2p)
legend('Signal R2 and Deconvolved');
ylabel('Power');
xlabel('Time(us)');
grid minor;
% Show specters
[output1, as_input1] = spect(R1, Fs, 'Plot', [true,false]);
[output2, as_input2] = spect(R2, Fs, 'Plot', [true,false]);
[Deconvolved1, as_input3] = spect(R1c, Fs, 'Plot', [true,false]);
[Deconvolved2, as_input4] = spect(R2c, Fs, 'Plot', [true,false]);
[chirp, as_input5] = spect(ch,Fs, 'Plot', [true,false]);
figure;
plot (output1/1e6, as_input1./max(as_input1(:)),output2/1e6, as_input2./max(as_input2(:)), Deconvolved1/1e6, as_input3./max(as_input3(:)),Deconvolved2/1e6, as_input4./max(as_input4(:)), chirp/1e6, as_input5./max(as_input5(:)));
title('Normalized frequency spectrum soon after reading signals');
legend('output1','output2','Deconvolved1','Deconvolved2','chirp');
set(gca, 'XLim', [0, 2]);
grid minor;
ylabel('Power');
xlabel('Frequency (MHz)');
%% THE DOUBT IS WHY THE SPECTRUM AFTER DECONVOLUTION IS SO DIFFERENT OF BEFORE.
Thanks for reading.
Tiago Dourado