Forums

anyone have an answer to this?

Started by bronx April 7, 2008
On Apr 8, 2:02&#2013266080;pm, maury <maury...@core.com> wrote:
> On Apr 7, 4:31&#2013266080;pm, "bronx" <branko_blagoje...@hotmail.com> wrote: > > > > > > > >You should post the equations for the frequencies, length of signals, > > >size of FFT, ... > > > >I suspect this method is going to run into problems where there are a > > >lot of FFT bins where there is virtually no signal, such that the > > >phase is really based more on numerical noise than any signal > > >present. &#2013266080;In that case the phase difference between these bins would > > >not be dependent on the delay and would not contribute to your > > >result. &#2013266080;I don't know if you have this problem or not, but it would > > >somethng to check for. > > > >Dirk > > > here is my code in detail > > > freq1 = 0.5e3; &#2013266080; &#2013266080; %Frequency of incoming source in Hz > > time = 20e-3; &#2013266080; %Total time of simulation in s > > tau1 = 0.375e-3; &#2013266080; &#2013266080; %time difference in arrival of 2 signals > > > N = pow2(11); > > %t = (0:N-1)*time; > > t = linspace(0,time,N); > > > % Signals arriving at input of microphones...with different phase > > > % MIC A > > x1 = sin(2*pi*freq1*t); > > > % MIC B > > x2 = sin(2*pi*freq1*(t + tau1)); > > > figure(1); > > plot(t, x1, "r"); > > hold on; > > plot(t, x2, "g"); > > > x1pad = [(x1) zeros(1,N-1)]; > > x2pad = [(x2) zeros(1,N-1)]; > > > NPC = length(x1pad); > > X1 = fft(x1pad); > > X2 = fft(x2pad); > > XR_CC_num = X1.*conj(X2); > > XR_CC_den = abs(XR_CC_num); > > XR = (XR_CC_num)./(XR_CC_den); > > %Phase Correlation > > PhaseCorr = fftshift((ifft(XR))); > > %Cross Correlation > > CrossCorr = fftshift(ifft(XR_CC_num)); > > CorrTime = time*linspace(-1,1, NPC); > > > figure(2); > > %Plot Phase correlation in ms > > plot(1e3*CorrTime,real(PhaseCorr),'r'); > > > figure(3); > > %Plot Cross Correlation in ms > > plot(1e3*CorrTime,real(CrossCorr),'g'); > > > %the location of the peak of each result is the time delay between signal > > %%x1 and signal x2 > > > let me know if you know why the peaks don't show the same time delay > > index. > > > Thanks > > x1 is a delayed version of x2 by tau1. &#2013266080;However, the question I think > you should be asking yourself is why you think the FFT of x1 would be > different than the FFT of x2. > > Maurice Givens- Hide quoted text - > > - Show quoted text -
Do you think they should be the same? You will notice that the OP never takes the FFT of x1 and x2. Dirk
bronx wrote:
>> You should post the equations for the frequencies, length of signals, >> size of FFT, ... >> >> I suspect this method is going to run into problems where there are a >> lot of FFT bins where there is virtually no signal, such that the >> phase is really based more on numerical noise than any signal >> present. In that case the phase difference between these bins would >> not be dependent on the delay and would not contribute to your >> result. I don't know if you have this problem or not, but it would >> somethng to check for. >> >> Dirk >> > > here is my code in detail > > freq1 = 0.5e3; %Frequency of incoming source in Hz
. . .
> > let me know if you know why the peaks don't show the same time delay > index. > > Thanks
bronx, I had seen this in a few years ago in a toy radar problem. So I thought I'd stab at this one. Two main points: noise and choice of signal. =Noise= At first I thought it might be due to some numerical issue as Dirk suggested so I added some noise, but to no avail. I left the noise code in but notice the magnitude is set to a very small level compare to the signal. Some noise is necessary to avoid any dividing by too small of numbers. Try it with zero and see the difference. Also experiment with larger noise to see how the accuracy decreases. =Signal selection= I believe the main issue with your example is the signal. From my understanding you need a somewhat richer signal than a perfect cosine. What's a simple signal with rich frequency content? A pulse, which is what some simple systems use (I've seen it in cheap ultrasound). Better than that would be a chirp signal. If you are looking at audio source location, most audio should be rich enough. =Misc= On a side note, for efficiency it would be ideal if you made the FFT lengths a power of two. Currently they are an odd length, although it looks you were aiming at that with the definition of N. Cheers, J Elms PS: Maybe I shouldn't just give you the answer, but I'm feeling generous. Modified code: freq1 = 0.5e3; %Frequency of incoming source in Hz time = 20e-3; %Total time of simulation in s tau1 = 0.375e-3; %time difference in arrival of 2 signals N = pow2(11); %t = (0:N-1)*time; t = linspace(0,time,N); % Signals arriving at input of microphones...with different phase NoiseMag = 1e-5; tPulseLen = 1e-3; % 1ms length pulse nPulseLen = round(tPulseLen/t(2)); % convert to samples BaseSig1 = ( (t>0) & (t<0+tPulseLen) ); BaseSig2 = ( (t>tau1) & (t<tau1+tPulseLen) ); % MIC A x1 = BaseSig1 + NoiseMag*rand(1,N); % MIC B x2 = BaseSig2 + NoiseMag*rand(1,N); %delay % figure(1); % plot(t, x1, "r"); % hold on; % plot(t, x2, "g"); x1pad = [(x1) zeros(1,N)]; x2pad = [(x2) zeros(1,N)]; NPC = length(x1pad); X1 = fft(x1pad); X2 = fft(x2pad); XR_CC_num = X1.*conj(X2); XR_CC_den = abs(XR_CC_num); XR = (XR_CC_num)./(XR_CC_den); %Phase Correlation PhaseCorr = real(fftshift((ifft(XR)))); %Cross Correlation CrossCorr = real(fftshift(ifft(XR_CC_num))); CorrTime = time*linspace(-1,1, NPC); % figure(2); % %Plot Phase correlation in ms % plot(1e3*CorrTime,real(PhaseCorr),'r'); % figure(3); % %Plot Cross Correlation in ms % plot(1e3*CorrTime,real(CrossCorr),'g'); %the location of the peak of each result is the time delay between signal %%x1 and signal x2 [tmp,indPhase] = max(PhaseCorr); [tmp,indCross] = max(CrossCorr); tPhase = CorrTime(indPhase) tCross = CorrTime(indCross) figure; %received and padded signals subplot(3,1,1) plot(x1pad ) hold on plot(x2pad ) legend('x1pad','x2pad') subplot(3,1,2) plot(abs(X1),':') hold on plot(abs(X2),'--') legend('|FFT(x1pad)|','|FFT(x2pad)|') subplot(3,1,3) plot(CrossCorr/max(CrossCorr) ) hold on plot(PhaseCorr/max(PhaseCorr) ) legend('Normalized CrossCorr','Normalized PhaseCorr')
On Apr 9, 4:23&#2013266080;pm, J Elms <jeffrey.e...@gatech.edu> wrote:
> bronx wrote: > >> You should post the equations for the frequencies, length of signals, > >> size of FFT, ... > > >> I suspect this method is going to run into problems where there are a > >> lot of FFT bins where there is virtually no signal, such that the > >> phase is really based more on numerical noise than any signal > >> present. &#2013266080;In that case the phase difference between these bins would > >> not be dependent on the delay and would not contribute to your > >> result. &#2013266080;I don't know if you have this problem or not, but it would > >> somethng to check for. > > >> Dirk > > > here is my code in detail > > > freq1 = 0.5e3; &#2013266080; &#2013266080; %Frequency of incoming source in Hz > > . > . > . > > > > > let me know if you know why the peaks don't show the same time delay > > index. > > > Thanks > > bronx, > > I had seen this in a few years ago in a toy radar problem. So I thought > I'd stab at this one. Two main points: noise and choice of signal. > > =Noise= > At first I thought it might be due to some numerical issue as Dirk > suggested so I added some noise, but to no avail. I left the noise code > in but notice the magnitude is set to a very small level compare to the > signal. Some noise is necessary to avoid any dividing by too small of > numbers. Try it with zero and see the difference. Also experiment with > larger noise to see how the accuracy decreases. > > =Signal selection= > I believe the main issue with your example is the signal. &#2013266080;From my > understanding you need a somewhat richer signal than a perfect cosine. > > What's a simple signal with rich frequency content? > A pulse, which is what some simple systems use (I've seen it in cheap > ultrasound). Better than that would be a chirp signal. If you are > looking at audio source location, most audio should be rich enough. > > =Misc= > On a side note, for efficiency it would be ideal if you made the FFT > lengths a power of two. &#2013266080;Currently they are an odd length, although it > looks you were aiming at that with the definition of N. > > Cheers, > J Elms > > PS: Maybe I shouldn't just give you the answer, but I'm feeling generous. > > Modified code: > > freq1 = 0.5e3; &#2013266080; &#2013266080; %Frequency of incoming source in Hz > time = 20e-3; &#2013266080; %Total time of simulation in s > tau1 = 0.375e-3; &#2013266080; &#2013266080; %time difference in arrival of 2 signals > > N = pow2(11); > %t = (0:N-1)*time; > t = linspace(0,time,N); > > % Signals arriving at input of microphones...with different phase > NoiseMag &#2013266080;= 1e-5; > tPulseLen = 1e-3; % 1ms length pulse > nPulseLen = round(tPulseLen/t(2)); % convert to samples > > BaseSig1 = ( (t>0) & (t<0+tPulseLen) ); > BaseSig2 = ( (t>tau1) & (t<tau1+tPulseLen) ); > > % MIC A > x1 = BaseSig1 + NoiseMag*rand(1,N); > > % MIC B > x2 = BaseSig2 + NoiseMag*rand(1,N); &#2013266080;%delay > > % figure(1); > % plot(t, x1, "r"); > % hold on; > % plot(t, x2, "g"); > > x1pad = [(x1) zeros(1,N)]; > x2pad = [(x2) zeros(1,N)]; > > NPC = length(x1pad); > X1 = fft(x1pad); > X2 = fft(x2pad); > XR_CC_num = X1.*conj(X2); > XR_CC_den = abs(XR_CC_num); > XR = (XR_CC_num)./(XR_CC_den); > %Phase Correlation > PhaseCorr = real(fftshift((ifft(XR)))); > %Cross Correlation > CrossCorr = real(fftshift(ifft(XR_CC_num))); > CorrTime = time*linspace(-1,1, NPC); > > % figure(2); > % %Plot Phase correlation in ms > % plot(1e3*CorrTime,real(PhaseCorr),'r'); > > % figure(3); > % %Plot Cross Correlation in ms > % plot(1e3*CorrTime,real(CrossCorr),'g'); > > %the location of the peak of each result is the time delay between signal > %%x1 and signal x2 > > [tmp,indPhase] = max(PhaseCorr); > [tmp,indCross] = max(CrossCorr); > tPhase = CorrTime(indPhase) > tCross = CorrTime(indCross) > > figure; > > %received and padded signals > subplot(3,1,1) > plot(x1pad ) > hold on > plot(x2pad ) > legend('x1pad','x2pad') > > subplot(3,1,2) > plot(abs(X1),':') > hold on > plot(abs(X2),'--') > legend('|FFT(x1pad)|','|FFT(x2pad)|') > > subplot(3,1,3) > plot(CrossCorr/max(CrossCorr) ) > hold on > plot(PhaseCorr/max(PhaseCorr) ) > legend('Normalized CrossCorr','Normalized PhaseCorr')- Hide quoted text - > > - Show quoted text -
Not only is the signal richer in frequency content, but over the interval they are defined, your delayed signal is really the original signal delayed, which was not true for the OP's signal. That makes the FFTs the same except for an approximately linear phase term which is important for the method to work. It is not quite linear phase due to the essentially random phase near zeros in the spectrum, but most of the phase values should be correct. Dirk
dbell wrote:
> On Apr 9, 4:23 pm, J Elms <jeffrey.e...@gatech.edu> wrote: >> bronx wrote: >>>> You should post the equations for the frequencies, length of signals, >>>> size of FFT, ... >>>> I suspect this method is going to run into problems where there are a >>>> lot of FFT bins where there is virtually no signal, such that the >>>> phase is really based more on numerical noise than any signal >>>> present. In that case the phase difference between these bins would >>>> not be dependent on the delay and would not contribute to your >>>> result. I don't know if you have this problem or not, but it would >>>> somethng to check for. >>>> Dirk >>> here is my code in detail >>> freq1 = 0.5e3; %Frequency of incoming source in Hz >> . >> . >> . >> >> >> >>> let me know if you know why the peaks don't show the same time delay >>> index. >>> Thanks >> bronx, >> >> I had seen this in a few years ago in a toy radar problem. So I thought >> I'd stab at this one. Two main points: noise and choice of signal. >> >> =Noise= >> At first I thought it might be due to some numerical issue as Dirk >> suggested so I added some noise, but to no avail. I left the noise code >> in but notice the magnitude is set to a very small level compare to the >> signal. Some noise is necessary to avoid any dividing by too small of >> numbers. Try it with zero and see the difference. Also experiment with >> larger noise to see how the accuracy decreases. >> >> =Signal selection= >> I believe the main issue with your example is the signal. From my >> understanding you need a somewhat richer signal than a perfect cosine. >> >> What's a simple signal with rich frequency content? >> A pulse, which is what some simple systems use (I've seen it in cheap >> ultrasound). Better than that would be a chirp signal. If you are >> looking at audio source location, most audio should be rich enough. >> >> =Misc= >> On a side note, for efficiency it would be ideal if you made the FFT >> lengths a power of two. Currently they are an odd length, although it >> looks you were aiming at that with the definition of N. >> >> Cheers, >> J Elms >> >> PS: Maybe I shouldn't just give you the answer, but I'm feeling generous. >> >> Modified code: >> >> freq1 = 0.5e3; %Frequency of incoming source in Hz >> time = 20e-3; %Total time of simulation in s >> tau1 = 0.375e-3; %time difference in arrival of 2 signals >> >> N = pow2(11); >> %t = (0:N-1)*time; >> t = linspace(0,time,N); >> >> % Signals arriving at input of microphones...with different phase >> NoiseMag = 1e-5; >> tPulseLen = 1e-3; % 1ms length pulse >> nPulseLen = round(tPulseLen/t(2)); % convert to samples >> >> BaseSig1 = ( (t>0) & (t<0+tPulseLen) ); >> BaseSig2 = ( (t>tau1) & (t<tau1+tPulseLen) ); >> >> % MIC A >> x1 = BaseSig1 + NoiseMag*rand(1,N); >> >> % MIC B >> x2 = BaseSig2 + NoiseMag*rand(1,N); %delay >> >> % figure(1); >> % plot(t, x1, "r"); >> % hold on; >> % plot(t, x2, "g"); >> >> x1pad = [(x1) zeros(1,N)]; >> x2pad = [(x2) zeros(1,N)]; >> >> NPC = length(x1pad); >> X1 = fft(x1pad); >> X2 = fft(x2pad); >> XR_CC_num = X1.*conj(X2); >> XR_CC_den = abs(XR_CC_num); >> XR = (XR_CC_num)./(XR_CC_den); >> %Phase Correlation >> PhaseCorr = real(fftshift((ifft(XR)))); >> %Cross Correlation >> CrossCorr = real(fftshift(ifft(XR_CC_num))); >> CorrTime = time*linspace(-1,1, NPC); >> >> % figure(2); >> % %Plot Phase correlation in ms >> % plot(1e3*CorrTime,real(PhaseCorr),'r'); >> >> % figure(3); >> % %Plot Cross Correlation in ms >> % plot(1e3*CorrTime,real(CrossCorr),'g'); >> >> %the location of the peak of each result is the time delay between signal >> %%x1 and signal x2 >> >> [tmp,indPhase] = max(PhaseCorr); >> [tmp,indCross] = max(CrossCorr); >> tPhase = CorrTime(indPhase) >> tCross = CorrTime(indCross) >> >> figure; >> >> %received and padded signals >> subplot(3,1,1) >> plot(x1pad ) >> hold on >> plot(x2pad ) >> legend('x1pad','x2pad') >> >> subplot(3,1,2) >> plot(abs(X1),':') >> hold on >> plot(abs(X2),'--') >> legend('|FFT(x1pad)|','|FFT(x2pad)|') >> >> subplot(3,1,3) >> plot(CrossCorr/max(CrossCorr) ) >> hold on >> plot(PhaseCorr/max(PhaseCorr) ) >> legend('Normalized CrossCorr','Normalized PhaseCorr')- Hide quoted text - >> >> - Show quoted text - > > Not only is the signal richer in frequency content, but over the > interval they are defined, your delayed signal is really the original > signal delayed, which was not true for the OP's signal. That makes the > FFTs the same except for an approximately linear phase term which is > important for the method to work. It is not quite linear phase due to > the essentially random phase near zeros in the spectrum, but most of > the phase values should be correct. > > Dirk
That's a good point, so I tried it with a real audio clip. For simulation purposes, I picked samples in an interlaced fashion so no samples were exactly overlapped between the signals (highly correlated, but not identical). I got accuracy within .25 ms and within 1 sample agreement between the two methods. I did encounter some anomalies with it detecting a zero delay when longer delays were used. But this begins to break the original assumption of this algorithm, which is that for the windows selected the data is highly correlated. However, even in these long delay cases there was also a significant peak at the correct location. An algorithm for peak detection could be used to improve results if these cases are encountered. Notes: With a real signal I did not have to add any noise. I also tried without padding the signals and with windowing the signal. These only had minor effects on the results. J Elms PS: Sorry if I ramble.
But why does a cross correlation plot contain information of the proper
delay between the signals, while the phase correlation does not?? If both
didn't give me anything, I would agree with you guys, but the fact that
one is giving something reasonable and the other is not makes me think
that there is something i am overlooking in the matlab coding, not in the
signals used.
>But why does a cross correlation plot contain information of the proper >delay between the signals, while the phase correlation does not?? If
both
>didn't give me anything, I would agree with you guys, but the fact that >one is giving something reasonable and the other is not makes me think >that there is something i am overlooking in the matlab coding, not in
the
>signals used. >
whoops, sorry guys, i didnt realize these responses came to 2 pages...this is actually a response to the last response on the first page....disregard it, and i will respond with more questions after i try your mod... thanks very much, bronx