Equation to Matlab code

12 ビュー (過去 30 日間)
Ajai Singh
Ajai Singh 2022 年 8 月 1 日
コメント済み: Ajai Singh 2022 年 8 月 2 日
Can someone please help to me understand how the equation (eq. 47) and code given in the picture relate ?
Any hint or suggestion would be of great help.
In the equation c is the order of the polynomial.
Thank you.
  5 件のコメント
Jeffrey Clark
Jeffrey Clark 2022 年 8 月 2 日
@Ajai Singh, @Torsten is correct - the LPC values are the computed multipliers for the x^(c-n) terms of the Legendre polynomials. So LPC(1,1) is the multiplier for the zero order Legendre polynomial, which is 1*x^0; similarly LPC(2,1:2) are the mulpliers for the 1st order polynomial 0*x^0 + 1*x^1. Check that the MATLAB code provided does compute the 2nd order polynomial multipliers when c = 3 to be LPC(3,1:3) = [-0.5 0 1.5] as in -0.5*x^0 + 0*x^1 + 1.5*x^2.
A table for the first 11 (0->10) can be seen in this article Legendre polynomials - Wikipedia
Note that MATLAB functions that operate on this short hand description of polynomials generally expects them in reverse order. See Polynomial roots - MATLAB roots (mathworks.com)
John D'Errico
John D'Errico 2022 年 8 月 2 日
編集済み: John D'Errico 2022 年 8 月 2 日
@jessupj - The Legendre polynomials from this family will be orthognal on the interval [-1,1], NOT [0,1]. In fact though, they could be evaluated on ANY interval, as they are still just polynomials. Some people choose to define the Legendre polynomials on the interval [0,1], which is entirely doable. But were you to do that, you would find a different family of polynomials. Regardless, you do not need to define what x is before you create any polynomials.

サインインしてコメントする。

採用された回答

John D'Errico
John D'Errico 2022 年 8 月 2 日
編集済み: John D'Errico 2022 年 8 月 2 日
(I really hope I am not doing your homework here, but the question is not posed as if it is, even though I assume you are a student.)
My take? Calling that a translation is a fairly liberal use of the word. I would also point out that it orders the coefficients of the polynomials in the wrong way. The MATLAB standard is to start with the highest order term FIRST.
On the interval [-1,1] the family of Legendre polynomials are an orthogonal family. We can find them listed in this link:
You can see the same 3 term recurrence relation posed in your question as Bonnet's recursion formula in that link.
The first few such polynomials are:
p_0(x) = 1
p_1(x) = x
p_2(x) = (3*x^2 - 1)/2
As you can see, the first two rows of LPC give us the first two polynomials, since we will interpret the coefficients in those rows as the coefficients of the powers of x. Essentially, those first two rows start the recursion.
But now, suppose I decided to write a truly direct translation of the recurrence relation? What would it look like? I might do this:
C = 7;
LPC = zeros(C+1,C+1);
LPC(1,1) = 1;
LPC(2,1:2) = [0 1];
for c = 1:C-1
LPC(c+2,1:c+2) = ((2*c+1)*[0,LPC(c+1,1:c+1)] - c*[LPC(c,1:c),0 0])/(c+1);
end
For example, do you see that multiplying the polynomial represented by coefficients [0 1], by x is represented by the code fragment where we shift the coefficients over by pre-pending a zero to the vector? So [0 1] to represent x would turn into [0 0 1], representing x^2.
LPC
LPC = 8×8
1.0000 0 0 0 0 0 0 0 0 1.0000 0 0 0 0 0 0 -0.5000 0 1.5000 0 0 0 0 0 0 -1.5000 0 2.5000 0 0 0 0 0.3750 0 -3.7500 0 4.3750 0 0 0 0 1.8750 0 -8.7500 0 7.8750 0 0 -0.3125 0 6.5625 0 -19.6875 0 14.4375 0 0 -2.1875 0 19.6875 0 -43.3125 0 26.8125
We can make that prettier and more easily read using syms.
syms x
LPC*x.^(0:C).'
ans = 
As you can see, this does generate the table of Legendre polynomials. We can see the same table in the wiki link I gave. And the code I wrote is a fairly faithful translation of the recurrence relation. It should be pretty easy to read, once you understand that multiplying a polynomial by x is equivalent to shifting the coefficients over and pre-pending a zero.
Perhaps a more interesting question is, why does the code you show also create the same set of coefficients? An analogy here would be, if I gave you the formula for a convolution of two functions, as a convolution. But suppose I wrote code using a Fourier transform, where I take the transforms of each function, multiply them in the Fourier domain, and then take the inverse Fourier transform? The result would still be the convolution of the functions, but I would surely not call it a literal translation of a convolution.
Anyway, why does the code written in your question work? You should see that I used only a single loop. Inside the loop, I extracted the coefficients of an entire lower order polynomial at once, then used them directly as a vector of coefficients. The loop you had is a double loop, but I'll claim it is a terribly poor implementation, even at that. Inside the inner loop, they first define LPC(i,j+1), and THEN define LPC(i,j), each time in the loop. UGH. Sorry, but that is just a convoluted code. What it does is allow you to not need a special case for the constant coefficient. I'll re-write my own code, now using a double, nested loop, DIRECTLY from the formula. As you should see, I create the constant term outside the loop.
C = 7;
LPC = zeros(C+1,C+1);
LPC(1,1) = 1;
LPC(2,1:2) = [0 1];
for c = 1:C-1
LPC(c+2,1) = -c*LPC(c,1)/(c+1); % the constant coefficient is assigned outside the loop
for j = 1:c+1
LPC(c+2,j+1) = ((2*c+1)*LPC(c+1,j) - c*LPC(c,j+1))/(c+1);
end
end
LPC*x.^(0:C).'
ans = 
As you should see, this does generate the same set of polynomials, using a nested loop.
  3 件のコメント
John D'Errico
John D'Errico 2022 年 8 月 2 日
The code you had was difficult to read, because while it was working, it did the job in a strange way. The simpler approach, where I deal with the complete polynomials in a single loop as vectors of polynomial coefficients is I think easier to follow. Easy to read code is always the best code, because then you can read it and debug it and most importantly, understand it.
Ajai Singh
Ajai Singh 2022 年 8 月 2 日
I agree , we referred to a papers supplimentary material and got that piece of code. The way i figured it out was comparing the code and the given equation and then solved c in terms of i and then everything was clear. Thank you for the suggestion. It made the code more readable. I appreciate it.

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeLinear Algebra についてさらに検索

製品


リリース

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by