clear all
Fs = 48000; % Sampling frequency
T = 1/Fs; % Sample time
L = 4096; % Length of signal
t = (0:L-1)*T; % Time vector
% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
a1 = 0.7;
a2 = 0.0;
f1 = 997;
f2 = 2700;
snr = -48; % db = 20 log (a/b)
b = 10^(snr/20)/((a1+a2)/2);
x1 = a1*sin(2*pi*f1*t) + a2*sin(2*pi*f2*t);
y1 = x1 + b*randn(size(t)); % Sinusoids plus noise
phase1 = .1;
del1 = 2*pi*phase1;
phase2 = .6;
del2 = 2*pi*phase2;
x2 = a1*sin(2*pi*f1*t+del1) + a2*sin(2*pi*f2*t+del2);
y2 = x2 + b*randn(size(t)); % Sinusoids plus noise
figure(1);
subplot(2,2,1);
plot(Fs*t(1:250),y1(1:250))
title('Signal Corrupted with Zero-Mean Random Noise y1(t)')
xlabel('time (milliseconds)')
subplot(2,2,3);
plot(Fs*t(1:250),y2(1:250))
title('Signal Corrupted with Zero-Mean Random Noise y2(t)')
xlabel('time (milliseconds)')
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y1 = fft(hamming(L)'.*y1,NFFT)/L;
Y2 = fft(hamming(L)'.*y2,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
% Plot single-sided amplitude spectrum.
subplot(2,2,2);
plot(f,2*abs(Y1(1:NFFT/2+1)))
title('Single-Sided Amplitude Spectrum of y1(t)')
xlabel('Frequency (Hz)')
ylabel('|Y1(f)|')
subplot(2,2,4);
plot(f,2*abs(Y2(1:NFFT/2+1)))
title('Single-Sided Amplitude Spectrum of y2(t)')
xlabel('Frequency (Hz)')
ylabel('|Y2(f)|')
alpha = 1;
rphat_nom = Y1.*conj(Y2);
%rphat_den = 1;
rphat_den = max(abs(rphat_nom), 1e-6);
%rphat_den = abs(rphat_nom);
rphat = rphat_nom./rphat_den;
TDOA = ifft(rphat);
maxTDOA = max(TDOA)
figure(2);
subplot(1,3,1);
plot(abs(TDOA(1:20)))
subplot(1,3,2);
plot(abs(TDOA(1:50)))
subplot(1,3,3);
plot(abs(TDOA));
评论1