Intégration d'une droite définie par spline.

Bonjour à tous,
Je viens ici vous exposer mon problème.

J'ai une série de point sur MATLAB correspondant à l'évolution du débit pendant un cycle cardiaque.
Le signal étant discret, je l'ai interpolé grâce à la fonction spline, jusqu'ici tout va bien.
Le problème est que j'ai besoin pour des calculs de post-traitement d'intégrer cette courbe, et je ne sais pas vraiment comment faire.

J'ai réalisé un petit script "test", celui-ci permet de dériver la courbe (le contraire de ce que je veux).

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13

clear all; close all; clc   X = 0:0.4:4*pi; Y = sin(X); S = spline(X,Y); M = diag(3:-1:1,1); S1 = S; S1.coefs = S1.coefs*M; x = linspace(X(1),X(end),1001); hold all plot(X,Y,'o-') plot(x,ppval(S,x),'r') plot(x,ppval(S1,x),'g')

---

Pour ce calcul j'utilise une matrice de transition M.

Ma question est, connaissez vous la matrice de transition permettant l'intégration ?
Ou bien, avez vous un moyen d'intégrer une courbe définie par spline ?

Je vous remercie de vos réponses

Ca me parait un peu compliqué, pourquoi ne fais tu pas une spline sur l'intégrale de tes valeurs?

Code:

---

1 2 3 4 5 6 7 8 9 10

X = 0:0.4:4*pi; Y = sin(X); S = spline(X,Y); Y_integ = cumsum(Y)*(X(2)-X(1)); S1 = spline(X,Y_integ); x = linspace(X(1),X(end),1001); hold all plot(X,Y,'o-') plot(x,ppval(S,x),'r') plot(x,ppval(S1,x),'g')

---

Le problème est que le calcul d'intégrale (ou plutôt de primitive) d'un signal discret n'est pas évident sous matlab.

Une chose pourrait m'aider, savais vous comment une courbe est reconstruite à partir des données en sorties de la fonction spline (breaks et coefs ?)

spline est une interpolation cubique.

Du coup tes coeffs S.coefs sont les coefficients d'un polynôme définit par morceaux de type aX^3 + bX^2+cX+d



C'est pourquoi ta matrice M fait bel et bien ta dérivation : fais bien

[aX^3 bX^2 cX d] * M = [0 3aX^2 2bX c];


le problème avec l'intégration, c'est que tu passes du degré 3 au degré 4, du coup il te faut un coeff de plus.
Tu peux t'en sortir en faisant un peut la même chose :

-tu rajoutes une colonne de 0 devant (et oui, avant d'intégrer, ton terme de degré 4 est nul)

Code:

---

1 2 3 4

S2 = S; S2.coefs = [zeros(size(S2.coefs,1),1) S2.coefs] S2.coefs = S2.coefs*M2;

---

en effet, ceci

[0 aX^3 bX^2 cX d] * M2 = [1/4aX^4 1/3bX^3 1/2cX^2 dX Constante];

avec

Code:

---

1 2 3 4 5 6 7 8 9 10

>> M2   M2 =   0 0 0 0 0 0.2500 0 0 0 0 0 0.3333 0 0 0 0 0 0.5000 0 0 0 0 0 1.0000 0

---

qui est donc une matrice 5*5 puisque tu dois passer au degré 4, réalise bien une primitive.

Il te reste à calculer la constante pour que ta primitive prenne la valeur voulue à l'origine (ici, elle s'annule en 0).

après t'aura quand même un problème avec ppval je pense.

J'ai essayé cela mais je n'ai aucune condition pour déterminer les constantes.

Savez vous comment on peut reconstruire la courbe spline à partir des coefficients calculés (pp.coefs) sans utiliser la fonction ppval, autrement dit à la "main" ?

Salut,

essaye ceci :

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35

X = 0:0.4:4*pi; Y = sin(X); S = spline(X,Y); M = diag(3:-1:1,1); M2 = zeros(5,5); M2(2,1) = 1/4; M2(3,2)= 1/3; M2(4,3) = 1/2; M2(5,4) = 1; S1 = S; S2 = S; S1.coefs = S1.coefs*M; S2.coefs = [zeros(size(S2.coefs,1),1) S2.coefs]; S2.coefs = S2.coefs*M2; % ajout des constantes (ceci est faux) S2.coefs(:,5) = pi^2/log(3212); x = linspace(X(1),X(end),1001); % on evalue la valeur des polynomes interpolateurs cubiques sur l'ensemble des points du domaine "serré" A = ppval(S,x); %on fait somme f(x)dt, proche de l'intégrale B = cumsum(A)*(x(2)-x(1)); %tracé fonction de départ hold all plot(X,Y,'o-') %tracé premiere spline plot(x,ppval(S,x),'r') %tracé dérivée plot(x,ppval(S1,x),'g') %tracé intégrale avec cumsum plot(x,B,'y') %tracé de l'intégrale calculée proprement de la première spline plot(x,x.^4*S2.coefs(1,1)+x.^3*S2.coefs(1,2)+x.^2*S2.coefs(1,3)+x.^1*S2.coefs(1,4),'m') ylim([-1.5 2.5])

---

L'idée de tes splines avec les breaks, c'est que tu approxime chaque portion de courbe entre deux "breaks" par un polynôme de degré 3. tu as 32 breaks, donc au final 31 polynomes de degré 3, dont les coefficients sont contenus dans S1.coefs.

au final, grossièrement, on peut dire que S1 représente une fonction définie par morceaux dans chaque intervalle [breakpoint breakpoint+1]. et chaque morceau est un polynome de degré 3.

ppval évalue la valeur de ta fonction S1 sur l'ensemble des poins du domaine x, soit 1001 points.

A partir de là tu as deux approches : soit tu fais l'intégration comme je te l'avait indiquée dans mon précédent post, et tu as des polynomes de degré 4 ou tu dois trouver la constante.

soit tu fais l'intégrale discrete avec la fonction cumsum dont je t'invite à aller lire la doc.


dans ce que je t'ai montré, j'ai tracé "à la main" le premier polynome correspondant à l'intégrale S2 de la spline S entre les premiers breakpoints (la courbe violette)

Code:

---

1 2 3

%tracé de l'intégrale calculée proprement de la première spline plot(x,x.^4*S2.coefs(1,1)+x.^3*S2.coefs(1,2)+x.^2*S2.coefs(1,3)+x.^1*S2.coefs(1,4),'m')

---

j'ai aussi tracé en jaune l'intégrale obtenue avec cumsum des 1001 points des valeurs de la spline (et pas que des points de départ comme je l'avais fait dans mon premier post).

Code:

---

1 2 3 4 5 6

% on évalue la valeur des polynômes interpolateurs cubiques sur l'ensemble des points du domaine "serré" A = ppval(S,x); %on fait somme f(x)dt, proche de l'intégrale B = cumsum(A)*(x(2)-x(1)); %tracé intégrale avec cumsum plot(x,B,'y')

---

tu constateras que sur l'intervalle [0 0.8] ou la courbe violette est valable, l'écart est très faible.

Personnellement, je pense que la deuxième solution est satisfaisante car précise et beaucoup plus simple à mettre en œuvre.

Effectivement ta solution est efficace et je t'en remercie.
Cependant comme je l'ai énoncé au début je fait des calculs à partir des débits lors d'un cycle cardiaque, je n'ai donc pas de conditions aux limites précise et valable à tout les coups.
Pour régler cela j'ai donc décidé de passer par une analyse de Fourier, sachant que l'intégration d'une série de fourier est très facile !

Je te remercie encore une fois de ton aide ;)