I had recently posted a question on applying powers on a sum of Gaussians (here) to enhance signal resolution artificially. The discussion has given rise to another query. Consider this problem. We have two Gaussian time functions G1 and G2. Suppose we wish sum these Gaussians and raise it to a power of 2.
Situation no. 1: (G1+G2)2, we should get another function G= G21 + 2G1G2 + G22. The product of two G1G2 is another Gaussian. We should see three peaks in this case.
Situation no. 2: We sum the two Gaussians element wise i.e., Si = g(t)1i + g(t)2i where g(t)i 's are the corresponding elements of G1 + G2. Now each summed element Si is raised to power 2. We get only two peaks as a result.
I discussed this issue with some of my colleagues including the person who had proposed this method; they seek agree that both operations are equivalent. It seems they are not, because in situation 2, we will never see G1G2. What would the discrete version of situation 1?
I feel there is some fallacy here and some mathematical insight should be able to settle this problem.
Thanks.
For Situation no. 2: Here is the MATLAB code.
t=[0:1/200:60]';
u1= 19; % mu the mean time peak i
u2= 22; % mu the mean time peak j
A= 250; %Area
s_ij=[0.84 0.87];; %s_ij represent the sigma of i,j pair; Choose elements by indexing
G_ij=A*normpdf(t, u1, s_ij(1,1))+ A*normpdf(t, u2, s_ij(1,2)); %Sum normpdf is a built in function for a unit area Gaussian from statistical package.
Power= (G_ij).^2;
subplot(2,1,1); plot(t,G_ij);title('Unresolved Gaussians');
subplot(2,1,2); plot (t, Power);title ('Gaussian Sum Raised to Power 2');
Answer
([a b c] .+ [d e f]).^2 = [a+d b+e c+f].^2
= [(a+d)^2 (b+e)^2 (c+f)^2]
= [a^2+2ad+d^2 b^2+2be+e^2 c^2+2cf+f^2]
I suspect that:
- The problem may be in your assumptions: "We should see three peaks in this case" sounds like a hypothesis that may be false.
- You may have a bug in your code, which might be why "We get only two peaks as a result".
Here is some code to illustrate that (A.+B).^2 = A.^2 .+ 2A.*B .+ B.^2
:
t=0:1/200:60;
gauss1 = 250*normpdf(t,19,0.84);
gauss2 = 250*normpdf(t,22,0.87);
m1 = gauss1.^2+2*gauss1.*gauss2+gauss2.^2;
m2 = (gauss1+gauss2).^2;
fprintf('Difference: %1.2f\n', sum(abs(m1-m2)));
subplot(2,1,1);plot(m1);title('Method 1');
subplot(2,1,2);plot(m2);title('Method 2');
I get:
Difference: 0.00
and
So, the math is correct. You may still be wondering, how come the third Gaussian G1G2 is not apparent in the plot? This is the point where it's useful to challenge your assumption that the third Gaussian should be visible. If we plot G1G2 by itself we obtain:
plot(2*gauss1.*gauss2);
Look at the amplitude and the position in the x-axis: this third Gaussian is not apparent in the plot because it is completely swamped by the other two, much larger, surrounding peaks.
No comments:
Post a Comment