Prev: Regarding Digital NCO design
Next: Alternatives to CEVA, Tensilica and CoolFlux NXP
From: Andor on 24 Sep 2009 10:46
On 24 Sep., 15:49, "sofiyya" <karimae...(a)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...(a)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...(a)gmail.com> wrote:
> >> >> On 24 Sep., 09:49, "sofiyya" <karimae...(a)gmail.com> wrote:
>
> >> >> > >On Sep 23, 7:25 am, "sofiyya" <karimae...(a)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) ...

From: sofiyya on 24 Sep 2009 11:12
>On 24 Sep., 15:49, "sofiyya" <karimae...(a)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...(a)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...(a)gmail.com> wrote:
>> >> >> On 24 Sep., 09:49, "sofiyya" <karimae...(a)gmail.com> wrote:
>>
>> >> >> > >On Sep 23, 7:25 am, "sofiyya" <karimae...(a)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) ...
>
>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...
From: sofiyya on 24 Sep 2009 11:20
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?
From: Dilip Warrier on 24 Sep 2009 15:30
On Sep 24, 3:49 am, "sofiyya" <karimae...(a)gmail.com> wrote:
> >On Sep 23, 7:25 am, "sofiyya" <karimae...(a)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_fft
array([ 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_fft
array([ 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)
>>> w
array([-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.

From: Andor on 24 Sep 2009 15:54
On 24 Sep., 15:49, Rune Allnor <all...(a)tele.ntnu.no> wrote:
> On 24 Sep, 15:42, Andor <andor.bari...(a)gmail.com> wrote:
>
>
>
>
>
> > On 24 Sep., 15:37, Rune Allnor <all...(a)tele.ntnu.no> wrote:
>
> > > On 24 Sep, 15:04, Andor <andor.bari...(a)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
First  |  Prev  |  Next  |  Last
Pages: 1 2 3 4 5
Prev: Regarding Digital NCO design
Next: Alternatives to CEVA, Tensilica and CoolFlux NXP