On 24 Sep., 15:37, Rune Allnor <all...@tele.ntnu.no> wrote:> On 24 Sep, 15:04, Andor <andor.bari...@gmail.com> wrote: > > > As I said, if you calculate 5 more terms in the w that I started > > below, I will tell you how to calculate w (given v1 and v2) in Matlab > > with one single command that has only 28 characters (challenge: who > > can do it in less?). > > I can do it with a 1-character command: > > %%%%%%%%%% File a.m %%%%%%%%%%%%% > > % implement computations here > > %%%%%%%%% End of a.m %%%%%%%%%%%% > > in which case the one-liner becomes > > >> a > > But that might not have been what you mean...?Nah, only Matlab and signal processing toolbox functions are allowed :-).
deconvolution problem
Started by ●September 23, 2009
Reply by ●September 24, 20092009-09-24
Reply by ●September 24, 20092009-09-24
On 24 Sep, 15:42, Andor <andor.bari...@gmail.com> wrote:> On 24 Sep., 15:37, Rune Allnor <all...@tele.ntnu.no> wrote: > > > > > > > On 24 Sep, 15:04, Andor <andor.bari...@gmail.com> wrote: > > > > As I said, if you calculate 5 more terms in the w that I started > > > below, I will tell you how to calculate w (given v1 and v2) in Matlab > > > with one single command that has only 28 characters (challenge: who > > > can do it in less?). > > > I can do it with a 1-character command: > > > %%%%%%%%%% File a.m %%%%%%%%%%%%% > > > % implement computations here > > > %%%%%%%%% End of a.m %%%%%%%%%%%% > > > in which case the one-liner becomes > > > >> a > > > But that might not have been what you mean...? > > Nah, only Matlab and signal processing toolbox functions are > allowed :-)Then you'll have to accept my solution above: I don't have the SPT. It would be much more interesting with only basic matlab commands allowed. Rune
Reply by ●September 24, 20092009-09-24
>As I said, if you calculate 5 more terms in the w that I started >below, I will tell you how to calculate w (given v1 and v2) in Matlab >with one single command that has only 28 characters (challenge: who >can do it in less?). > >:-) > > >On 24 Sep., 14:27, "sofiyya" <karimae...@gmail.com> wrote: >> Thank you for your response, >> >> Yes, with pencil and paper we can calculate w but it will be difficultfor>> a long vector, that's why I look for a methods to do it with matlab..Maybe>> it's impossible or maybe this depends on the coefficients of vect1 or >> vect2.. I guess that we can't always find w / conv(w,vect1)=vect2; >> >> >> >> >> >> >On 24 Sep., 10:45, Andor <andor.bari...@gmail.com> wrote: >> >> On 24 Sep., 09:49, "sofiyya" <karimae...@gmail.com> wrote: >> >> >> > >On Sep 23, 7:25 am, "sofiyya" <karimae...@gmail.com> wrote: >> >> >> > This is the result I get when I don't truncate w: >> >> >> > ans =3D >> >> >> > =A0 Columns 1 through 9 >> >> >> > =A0 =A0-0.0903 =A0 -0.2708 =A0 -0.5417 =A0 -0.9028 =A0 -0.7292=A0>> -0.5= >> >208 =A0 -0.2778 =A0 >> >> > -0.0000 =A0 =A00.2778 >> >> >> > =A0 Columns 10 through 16 >> >> >> > =A0 =A0 1.4583 =A0 =A02.7292 =A0 =A04.0903 =A0 =A05.5417 =A0 >> =A00.8333 = >> >=A0 =A00.5903 =A0 =A00.3125 >> >> >> > which is also wrong! >> >> >> Your method is flawed alltogehter. In general, w will have infinite >> >> length. Any truncation to use the FFT won't produce the expected >> >> results. w can be found through long division of v2 by v1. This is >> >> very simple to do with pencil and paper, the first couple of valuesin>> >> w are >> >> >> w =3D [1 0 0 0 -5 4 ...] >> >> >If you add 5 more terms yourself, I'll reveal the matlab one-liner >> >that calculates w to arbitrary length :-).- Zitierten Text ausblenden->> >> - Zitierten Text anzeigen - > >w must have only 6 items to get the right answer!
Reply by ●September 24, 20092009-09-24
>On 24 Sep, 15:04, Andor <andor.bari...@gmail.com> wrote: >> As I said, if you calculate 5 more terms in the w that I started >> below, I will tell you how to calculate w (given v1 and v2) in Matlab >> with one single command that has only 28 characters (challenge: who >> can do it in less?). > >I can do it with a 1-character command: > >%%%%%%%%%% File a.m %%%%%%%%%%%%% > >% implement computations here > >%%%%%%%%% End of a.m %%%%%%%%%%%% > >in which case the one-liner becomes > >>> a > >But that might not have been what you mean...? > >Rune >I know your function: is deconv(vect1,vect2);
Reply by ●September 24, 20092009-09-24
On 24 Sep., 15:49, "sofiyya" <karimae...@gmail.com> wrote:> >As I said, if you calculate 5 more terms in the w that I started > >below, I will tell you how to calculate w (given v1 and v2) in Matlab > >with one single command that has only 28 characters (challenge: who > >can do it in less?). > > >:-) > > >On 24 Sep., 14:27, "sofiyya" <karimae...@gmail.com> wrote: > >> Thank you for your response, > > >> Yes, with pencil and paper we can calculate w but it will be difficult > for > >> a long vector, that's why I look for a methods to do it with matlab.. > Maybe > >> it's impossible or maybe this depends on the coefficients of vect1 or > >> vect2.. I guess that we can't always find w / conv(w,vect1)=vect2; > > >> >On 24 Sep., 10:45, Andor <andor.bari...@gmail.com> wrote: > >> >> On 24 Sep., 09:49, "sofiyya" <karimae...@gmail.com> wrote: > > >> >> > >On Sep 23, 7:25 am, "sofiyya" <karimae...@gmail.com> wrote: > > >> >> > This is the result I get when I don't truncate w: > > >> >> > ans =3D > > >> >> > =A0 Columns 1 through 9 > > >> >> > =A0 =A0-0.0903 =A0 -0.2708 =A0 -0.5417 =A0 -0.9028 =A0 -0.7292 > =A0 > >> -0.5= > >> >208 =A0 -0.2778 =A0 > >> >> > -0.0000 =A0 =A00.2778 > > >> >> > =A0 Columns 10 through 16 > > >> >> > =A0 =A0 1.4583 =A0 =A02.7292 =A0 =A04.0903 =A0 =A05.5417 =A0 > >> =A00.8333 = > >> >=A0 =A00.5903 =A0 =A00.3125 > > >> >> > which is also wrong! > > >> >> Your method is flawed alltogehter. In general, w will have infinite > >> >> length. Any truncation to use the FFT won't produce the expected > >> >> results. w can be found through long division of v2 by v1. This is > >> >> very simple to do with pencil and paper, the first couple of values > in > >> >> w are > > >> >> w =3D [1 0 0 0 -5 4 ...] > > >> >If you add 5 more terms yourself, I'll reveal the matlab one-liner > >> >that calculates w to arbitrary length :-).- Zitierten Text ausblenden > - > > >> - Zitierten Text anzeigen - > > w must have only 6 items to get the right answer!You thinke w = deconv(vect1,vect2) gives you the correct answer? Try conv(deconv(vect1,vect2),vect1) ...
Reply by ●September 24, 20092009-09-24
>On 24 Sep., 15:49, "sofiyya" <karimae...@gmail.com> wrote: >> >As I said, if you calculate 5 more terms in the w that I started >> >below, I will tell you how to calculate w (given v1 and v2) in Matlab >> >with one single command that has only 28 characters (challenge: who >> >can do it in less?). >> >> >:-) >> >> >On 24 Sep., 14:27, "sofiyya" <karimae...@gmail.com> wrote: >> >> Thank you for your response, >> >> >> Yes, with pencil and paper we can calculate w but it will bedifficult>> for >> >> a long vector, that's why I look for a methods to do it withmatlab..>> Maybe >> >> it's impossible or maybe this depends on the coefficients of vect1or>> >> vect2.. I guess that we can't always find w / conv(w,vect1)=vect2; >> >> >> >On 24 Sep., 10:45, Andor <andor.bari...@gmail.com> wrote: >> >> >> On 24 Sep., 09:49, "sofiyya" <karimae...@gmail.com> wrote: >> >> >> >> > >On Sep 23, 7:25 am, "sofiyya" <karimae...@gmail.com> wrote: >> >> >> >> > This is the result I get when I don't truncate w: >> >> >> >> > ans =3D >> >> >> >> > =A0 Columns 1 through 9 >> >> >> >> > =A0 =A0-0.0903 =A0 -0.2708 =A0 -0.5417 =A0 -0.9028 =A0 -0.7292 >> =A0 >> >> -0.5= >> >> >208 =A0 -0.2778 =A0 >> >> >> > -0.0000 =A0 =A00.2778 >> >> >> >> > =A0 Columns 10 through 16 >> >> >> >> > =A0 =A0 1.4583 =A0 =A02.7292 =A0 =A04.0903 =A0 =A05.5417 =A0 >> >> =A00.8333 = >> >> >=A0 =A00.5903 =A0 =A00.3125 >> >> >> >> > which is also wrong! >> >> >> >> Your method is flawed alltogehter. In general, w will haveinfinite>> >> >> length. Any truncation to use the FFT won't produce the expected >> >> >> results. w can be found through long division of v2 by v1. Thisis>> >> >> very simple to do with pencil and paper, the first couple ofvalues>> in >> >> >> w are >> >> >> >> w =3D [1 0 0 0 -5 4 ...] >> >> >> >If you add 5 more terms yourself, I'll reveal the matlab one-liner >> >> >that calculates w to arbitrary length :-).- Zitierten Textausblenden>> - >> >> >> - Zitierten Text anzeigen - >> >> w must have only 6 items to get the right answer! > >You thinke w = deconv(vect1,vect2) gives you the correct answer? > >Try conv(deconv(vect1,vect2),vect1) ... > >No w = deconv(vect1,vect2) did not give the correct answer (give just afew correct items), that's why I said that maybe it's impossible...
Reply by ●September 24, 20092009-09-24
No w = deconv(vect1,vect2) did not give the correct answer (give just a few correct items), that's why I said that maybe it's impossible... You have another idea?
Reply by ●September 24, 20092009-09-24
On Sep 24, 3:49 am, "sofiyya" <karimae...@gmail.com> wrote:> >On Sep 23, 7:25 am, "sofiyya" <karimae...@gmail.com> wrote: > > This is the result I get when I don't truncate w: > > ans = > > Columns 1 through 9 > > -0.0903 -0.2708 -0.5417 -0.9028 -0.7292 -0.5208 -0.2778 > -0.0000 0.2778 > > Columns 10 through 16 > > 1.4583 2.7292 4.0903 5.5417 0.8333 0.5903 0.3125 > > which is also wrong!I think the difference is due to the difference between linear convolution and circular convolution. For discrete time sequences of finite range, the conv function is probably doing a linear convolution. However, the result that: Y[n]= X[n]H[n] only holds in the frequency domain if y[k] = circular_convolution(x[k],h[k]). My code below illustrates the point (it's in python, not matlab, but fairly easy to follow).>>> vect1=[1,2,3,4,5,6,7,8,9] >>> vect2=[1,2,3,4,0,0,0,0,0,0,0,0,0,0] >>> vect1_fft=pylab.fft(vect1,16) >>> vect1_fftarray([ 45.00000000 +0.j , -17.13707118-25.13669746j, 5.00000000 +9.65685425j, -5.61991440 -7.48302881j, 5.00000000 +4.j , -4.72323135 -3.34089319j, 5.00000000 +1.65685425j, -4.51978306 -0.99456184j, 5.00000000 +0.j , -4.51978306 +0.99456184j, 5.00000000 -1.65685425j, -4.72323135 +3.34089319j, 5.00000000 -4.j , -5.61991440 +7.48302881j, 5.00000000 -9.65685425j, -17.13707118+25.13669746j])>>> vect2_fft=pylab.fft(vect2,16) >>> vect2_fftarray([ 10.00000000+0.j , 6.49981314-6.58220534j, -0.41421356-7.24264069j, -4.05147161-2.43834568j, -2.00000000+2.j , 1.80883092+1.80429501j, 2.41421356-1.24264069j, -0.25717245-2.33956465j, -2.00000000+0.j , -0.25717245+2.33956465j, 2.41421356+1.24264069j, 1.80883092-1.80429501j, -2.00000000-2.j , -4.05147161+2.43834568j, -0.41421356+7.24264069j, 6.49981314+6.58220534j])>>> w_fft=[vect2_fft[i]/vect1_fft[i] for i in range(16)] >>> w_fft[(0.22222222222222221+0j), (0.058417319893240283+0.29840494835137932j), (-0.60895771340131999-0.27240496093971622j), (0.46832072266248786-0.18970249468587133j), (-0.04878048780487805+0.43902439024390244j), (-0.43535327530506751-0.074065018168567823j), (0.36086262044185607-0.36810749065626669j), (0.16291305636344114+0.48177921624409881j), (-0.40000000000000002+0j), (0.16291305636343922-0.48177921624410136j), (0.36086262044185735+0.3681074906562658j), (-0.43535327530506812+0.074065018168570293j), (-0.04878048780487805-0.43902439024390244j), (0.46832072266248881+0.18970249468587194j), (-0.60895771340132054+0.27240496093971661j), (0.058417319893240185-0.29840494835137948j)]>>> w=pylab.ifft(w_fft,16) >>> warray([-0.01643333 -8.67361738e-18j, -0.02079546 +1.46247015e-16j, 0.02610188 -2.66631703e-16j, 0.02402530 +7.92890601e-17j, 0.02227031 +1.38777878e-17j, 0.01976020 -4.95361808e-17j, 0.01672933 +9.89262816e-18j, -0.08658693 -1.90311363e-16j, -0.08000779 -5.37764278e-17j, -0.06939682 -1.20033039e-16j, -0.06005461 +4.75829945e-16j, 0.44817788 +7.92890601e-17j, 0.00533612 +4.85722573e-17j, 0.00647544 +2.33222052e-17j, -0.00283080 -2.19090871e-16j, -0.01054850 +3.17332423e-17j])>>> pylab.convolve(w,vect1)array([ -1.64333308e-02 -8.67361738e-18j, -5.36621200e-02 +1.28899780e-16j, -6.47890336e-02 -1.58525533e-19j, -5.18906475e-02 -4.99277708e-17j, -1.67219479e-02 -8.58192282e-17j, 3.82069538e-02 -1.71246867e-16j, 1.09865188e-01 -2.46781877e-16j, 9.49364871e-02 -5.12628249e-16j, -2.60208521e-16 -8.32251050e-16j, -3.67761377e-16 -1.18517072e-15j, -7.80625564e-16 -2.60279314e-15j, -2.91433544e-16 +4.14136563e-17j, -9.15933995e-16 -4.58383216e-16j, -1.27155231e-15 -3.60034221e-16j, -2.06605566e-15 +1.39485802e-16j, -2.58473798e-16 +1.25987159e-16j, 1.01643333e+00 +2.10463580e-15j, 2.05366212e+00 +2.90824645e-15j, 3.06478903e+00 +4.42819964e-15j, 4.05189065e+00 +1.09556024e-16j, 1.67219479e-02 -7.19508683e-16j, -3.82069538e-02 -1.32069442e-15j, -1.09865188e-01 -1.71795190e-15j, -9.49364871e-02 +2.85599181e-16j]) Note that this is linear convolution and hence is of size 16 + 9 -1 = 24. The relevant values of 1, 2, 3, and 4 show up at index 16, 17, 18, 19 resp. To get the circular convolution, do the following.>>> result=list(pylab.convolve(pylab.real(w),vect1)) >>> result[-0.016433330754891194, -0.053662120030136096, -0.064789033596004547, -0.051890647462305581, -0.016721947873467282, 0.038206953835679201, 0.10986518759613308, 0.094936487127274621, -2.6020852139652106e-16, -3.677613769070831e-16, -7.8062556418956319e-16, -2.9143354396410359e-16, -9.1593399531575415e-16, -1.2715523078909996e-15, -2.0660556598883772e-15, -2.5847379792054426e-16, 1.0164333307548912, 2.0536621200301362, 3.0647890335960049, 4.0518906474623053, 0.016721947873467952, -0.03820695383567816, -0.10986518759613122, -0.094936487127274566]>>> [result[i] + result[16+i] for i in range(16) if i < 8][1.0, 2.0, 3.0000000000000004, 3.9999999999999996, 6.6960326172704754e-16, 1.0408340855860843e-15, 1.8596235662471372e-15, 5.5511151231257827e-17] which gives you what you were looking for. The derivation of the circular convolution result from linear convolution is mentioned in most DSP textbooks. Dilip.
Reply by ●September 24, 20092009-09-24
On 24 Sep., 15:49, Rune Allnor <all...@tele.ntnu.no> wrote:> On 24 Sep, 15:42, Andor <andor.bari...@gmail.com> wrote: > > > > > > > On 24 Sep., 15:37, Rune Allnor <all...@tele.ntnu.no> wrote: > > > > On 24 Sep, 15:04, Andor <andor.bari...@gmail.com> wrote: > > > > > As I said, if you calculate 5 more terms in the w that I started > > > > below, I will tell you how to calculate w (given v1 and v2) in Matlab > > > > with one single command that has only 28 characters (challenge: who > > > > can do it in less?). > > > > I can do it with a 1-character command: > > > > %%%%%%%%%% File a.m %%%%%%%%%%%%% > > > > % implement computations here > > > > %%%%%%%%% End of a.m %%%%%%%%%%%% > > > > in which case the one-liner becomes > > > > >> a > > > > But that might not have been what you mean...? > > > Nah, only Matlab and signal processing toolbox functions are > > allowed :-) > > Then you'll have to accept my solution above: > I don't have the SPT. It would be much more > interesting with only basic matlab commands allowed.Interesting perhaps, but it is much more instructive to use the "filter" command. Consider the problem: given two vectors x and y, find w such that w*x = y. The solution is long division of y by x. Now remember that the impulse response of a rational transfer function filter is the long division of the numerator by the denominator polynomial. So w is simply the impulse response of the filter with numerator y and denominator x! This will, in general, be an IIR filter, thus w has infinite length. The neat Matlab one-liner thus is w = filter(y,x,[1 zeros(1,n)]) where n is the number of terms you want to compute in w. Regards, Andor
Reply by ●September 24, 20092009-09-24
On 24 Sep, 21:54, Andor <andor.bari...@gmail.com> wrote:> On 24 Sep., 15:49, Rune Allnor <all...@tele.ntnu.no> wrote: > > > > > > > On 24 Sep, 15:42, Andor <andor.bari...@gmail.com> wrote: > > > > On 24 Sep., 15:37, Rune Allnor <all...@tele.ntnu.no> wrote: > > > > > On 24 Sep, 15:04, Andor <andor.bari...@gmail.com> wrote: > > > > > > As I said, if you calculate 5 more terms in the w that I started > > > > > below, I will tell you how to calculate w (given v1 and v2) in Matlab > > > > > with one single command that has only 28 characters (challenge: who > > > > > can do it in less?). > > > > > I can do it with a 1-character command: > > > > > %%%%%%%%%% File a.m %%%%%%%%%%%%% > > > > > % implement computations here > > > > > %%%%%%%%% End of a.m %%%%%%%%%%%% > > > > > in which case the one-liner becomes > > > > > >> a > > > > > But that might not have been what you mean...? > > > > Nah, only Matlab and signal processing toolbox functions are > > > allowed :-) > > > Then you'll have to accept my solution above: > > I don't have the SPT. It would be much more > > interesting with only basic matlab commands allowed. > > Interesting perhaps, but it is much more instructive to use the > "filter" command. > > Consider the problem: given two vectors x and y, find w such that w*x > = y. The solution is long division of y by x. Now remember that the > impulse response of a rational transfer function filter is the long > division of the numerator by the denominator polynomial. So w is > simply the impulse response of the filter with numerator y and > denominator x! This will, in general, be an IIR filter, thus w has > infinite length. The neat Matlab one-liner thus is > > w = filter(y,x,[1 zeros(1,n)]) > > where n is the number of terms you want to compute in w.Why do you need the SPT for that one? Apart from that, your approach is sensitive to the usual questions about poles and stability. If the denomiator polynomial x has poles on or outside the unit circle your approach blows up. Rune
