How determine the delay in my signal practically
The Least Mean Square solution to find the "channel" or response of the filter is provided by the following MATLAB/Octave Code using  the input to the filter as tx and the output of the filter as rx. For more details on how this works, see this post: Compensating Loudspeaker frequency response in an audio signal:

function coeff = channel(tx,rx,ntaps)
    % Determines channel coefficients using the Wiener-Hopf equations (LMS Solution)
    % TX = Transmitted (channel input) waveform, row vector, length must be &gt;&gt; ntaps 
    % RX = Received (ch output) waveform, row vector, length must be &gt;&gt; ntaps 
    % NTAPS = Number of taps for channel coefficients
    % Dan Boschen 1/13/2020

    tx= tx(:)';   % force row vector
    rx= rx(:)';   % force row vector
    depth = min(length(rx),length(tx));
    A=convmtx(rx(1:depth).',ntaps);
    R=A'*A;       % autocorrelation matrix
    X=[tx(1:depth) zeros(1,ntaps-1)].';
    ro=A'*X;      % cross correlation vector
    coeff=(inv(R)*ro);   %solution
end


The case the OP uses of a comb filter is one of the most challenging as it depends on signal energy at each frequency for a solution (this is specifically why linear equalizers, which this function is doing if you swap rx and tx, do not perform well in frequency selective channels and result in noise enhancement at the null locations).  Below the frequency response of the test filter as determined with MATLAB or Octave showing the multiple frequency nulls associated with such a comb filter:

b=[1,zeros(1,39),-1];
freqz(b,1,2^14)     % 2^14 samples to show nulls




MATLAB or Octave script to demonstrate use of above function and determine the delay between output and input:

%% Filter with OP's example
b=[1,zeros(1,39),-1];     % numerator coefficients
a = 1;                    % denominator coefficients

%% Generate signal using OP's code
Fs=1000;
t=0:1/Fs:(5-1/Fs);
wi=blackman(length(t))';
rn=+rand(1,length(t))*.2;
x1=sin(2*pi*13*t).*wi +rn;
x2=sin(2*pi*25*t).*wi +rn;
x3=sin(2*pi*2*t).*wi +rn;
x=[x1,x2,x3];

% Filter
y=filter(b,a,x);

%% Test filter estimation
cf=channel(x,y,61);

%compare original and estimated channel
subplot(2,1,1)
stem(b)
title("Actual Channel Response")
xlabel("Sample Number")
subplot(2,1,2)

stem(cf)
title("Estimated Channel Response")
xlabel("Sample Number")




We could have used exactly 41 taps when calling the function channel and it would have resolved the solution properly; however it is best practice to start with a much longer filter length, evaluate the result and then reduce the taps accordingly. In actual practice under noise conditions using more taps than necessary will lead to noise enhancement, so it is good to do the final solution with the minimum taps necessary to capture the dominant tap weights. 

Observe with the groupdelay plot using MATLAB and Octaves grpdelay command the issue with not being able to resolve the delay where no signal gets through the filter (it would be difficult for anything to determine the delay of a single tone at one of these frequencies that is nulled by the filter!) but it is able to accurately determined the delay where signal energy exists. Similarly the waveform itself must have energy at all frequencies where we are seeking a solution. The spectral density of the OP's test waveform was sufficiently spread out over all frequencies to be suitable for this purpose (this is the reason why pseudo-random waveform make for good "channel sounding" patterns).

This plot is for comparison to the OP's plot showing that the group delay for this filter is 20 samples. 



Here is another test case showing how this works with a more reasonable channel (no deep nulls) with the following coefficients and frequency response below that:

b = [0.2 .4 -.3 .4 .3 .1];




The solution is indistinguishable between actual and estimate, so I added noise to make it more interesting, using same x and y from code above:

% add noise
noise = 0.351*randn(1,length(y));
yn = y + noise;
snr = 20*log10(std(yn)./std(noise));
%% Test filter estimation
cf=channel(x,yn,10);