Detrending algorithme et nombre de points

**darkwhite** · 10/08/2013, 04h16

Bonjour à tous,
J'essaie de réaliser un algo inspiré du passage de l'article suivant :

We have devised a new local adaptive detrending algorithm with these properties. It first partitions a time series into segments (or windows) of length 2n + 1 points, where neighboring segments overlap by n+1 points. For each segment, we fit a best polynomial of order K.Denote the fitted
polynomials for the ith and (i+1)th segments by y(i)(l1), y(i+1)(l2), l1, l2 = 1, . . . , 2n+1,
respectively. Note that the length of the last segment may be smaller than 2n + 1. We define the trend for the overlapped region as

y(c)(l)=(1-(l-1)/n)*y(i)(l+n) + (l-1)/n*y(i+1)(l), l=1,2,..,n+1

En fait on recouvre le signal par des segments de taille M=2n+1. Et chaque segments recouvre son voisin sur n+1 points.
Ensuite on crée le polynome Yc sur chaque intersection.
Ce qui m’étonne (voir m’énerve) c'est qu'au final le polynôme Yc que j'obtiens est plus long que le signal de base... Alors que normalement, avec la division euclidienne on devrait tomber sur une taille au moins plus petite.

Voici mon code :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
load 'input.mat'
n=50;
M=2*n+1;
 
Nseg=floor((length(input)-(2*n+1))/n);
 
K=2;
x=1:M;%taille d'un segment y
l=1:n+1;%taille d'un segments de recouvrement Yc
 
for i=1:Nseg
    seg(i,:)=input(1+(i-1)*n:(i-1)*n+M);
    P=polyfit(x,seg(i,:),K);
    y(i,:)=polyval(P,x);
end
 
%On construit 'the trend for the overlapped region' 
%Il s'applique sur chaque intersections de deux segments
 
for i=1:Nseg-1
Yc(i,l)=y(i,l+n).*(1-(l-1)/n)+y(i+1,l).*(l-1)/n;
end

seg : les segments qui recouvrent le signal input
y : estimation polynomiale sur chaque segment
Yc: estimation polynomiale sur les intersections de chaque segments

cependant la taille de Yc est plus grande que input ...

En gambergeant j'ai essayé de recenser
les points en trop:
*Le dernier point du ieme segment correspond au premier point du i+2e segment si on est d'accord.
les points qui manquent:
*Les L points constitués du nombre de points entre la fin du dernier segments et le nombre de point du signal.

Cependant ça ne fait pas le compte.

Voilà si vous avez la patience pour vous prendre la tête sur un problème de dénombrement comme celui là je vous serais extrêmement reconnaissant, sinon bonsoir !.
ps: je n'ai pas trouvé de balise latex, j'espere que c'est pas trop crade.

**tachmou** · 14/08/2013, 13h42

Bonjour,

Je n'ai pas tout bien lu dans le détail, mais que contient "input" stp ? Peux-tu faire un size(input) et donner la valeur ici stp. J'ai l'impression que tu travailles sur des tableaux alors qu'à la base tu utilises "length" qui te renvoies la taille d'une liste et non d'un tableau...

Par ailleurs, es-tu sûr de tes définitions de "x" et de "l" dans ton programme ? Tu dis en commentaire qu'ils contiennent des tailles alors que pour moi ce sont des listes allant de 1 à la taille dont tu parles...

Quoiqu'il en soit, avant de travailler sur des tableaux de dimension mxn, je te conseilles de travailler sur des listes de tailles 1xn...

Ce n'est pas du chipotage : j'essaye de comprendre ton programme... et comme je ne suis pas du domaine, bah je pige rien :/ En effet, sur les lignes dont je parle, je ne sais pas ce qui est juste : tableau ou liste ? Liste ou valeur ? Du coup, ça donne l'impression que tu te contredis sur la même ligne haha

Essaye de proposer un code complet en corrigeant ces "petits" détails... ça nous permettra de comprendre plus facilement, et peut-être que ça te permettra de te rendre compte de certaines erreurs par toi-même

Cordialement,

**darkwhite** · 14/08/2013, 16h43

Bonjour,
En tout cas merci de t'interesser à mon problème ...
Input est un signal colonne de 800 points :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
size(input)
ans =
   800     1

x et l sont inspiré du texte en anglais. x est l'index de tout polynome y
car on a dans l'article y(l1) avec l1=1,2,..,2*n+1=M ici l1=x. De même l correspond à l'index de tout polynome Yc car on a Yc(l) avec l=1,2,3,..,n.
Effectivement, le terme de taille est mal choisit.

Quoiqu'il en soit, avant de travailler sur des tableaux de dimension mxn, je te conseilles de travailler sur des listes de tailles 1xn...

A quel niveau ? pour stocker les y ? je ne vois pas vraimment comment faire avec un seul vecteur. La base du problème est de decouper un signal en signal de 2n+1 points. Mais chaque segment recouvre son voisin à moitié (n+1 points). J'ai mis un schema rapide pour décrire.

http://hpics.li/be67722
Voilà si ya besoin de plus de précisions j'en donne plus !

**tachmou** · 14/08/2013, 17h26

Okay, merci, j'avais un doute. En ce qui concerne ma dernière remarque, je pensais justement que ton tableau était de taille mxn (genre tu voulais traiter plusieurs signaux en parallèle)... donc ton length renvoyait totalement autre chose et que par la suite ça devenait n'importe quoi... sachant que je pensais que pour toi la taille se stockait dans une liste : je n'ai pas cherché plus loin. Bref j'étais à l'ouest : comme quoi ces "détails" peuvent mener loin ^^

Peux-tu donc poster un code en corrigeant ces détails et en commentant un poil plus (en étant rigoureux donc). Juste cette étape permettra à d'autres de comprendre plus facilement. Y compris moi

Conseil perso : x, tu peux le renommer index_segment ou i_segment ou i_seg et l tu peux le renommer index_recouvrement... ça parle plus. Cela alourdit le code mais peut le rendre tellement plus compréhensible. A titre informatif, j'ai développé des codes de près de 10 000 lignes... dont 8 000 de commentaires

==> soit 2 000 seulement, mais tellement efficaces et faciles à lire/relire/corriger, que ce soit immédiatement ou dans plusieurs années quand tu auras tout oublié ^^ Je pense notamment à certaines petites fonctions qui font à peine 5 lignes de code mais carrément près d'une centaine de commentaire

Bref, tu n'es pas obligé de suivre tous mes conseils hein mais reposte au moins un code concis, clair et sans ambiguïté.

**darkwhite** · 14/08/2013, 20h05

Merci, je comprends tout à fait les raison des commentaires.
Voici le code plus commenté:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
36
37
38
39
40
41
42
43
44
45
 
clc;
close all;
clear all;
 
 
load 'input.mat'
 
n=4;
%Taille d'un segment
M=2*n+1;
 
%Nombre de segments de taille M qui recouvrent le signal
%et qui ont n+1 points en commun
Nseg=floor(length(input)/n)-1;
 
 
K=2;%Degre du polynome y
index_y=1:M;%axe des abscisse de chaque polynome y
index_yc=1:n+1;%axe de abscisse de chaque Yc
 
 
B=[];
for i=1:Nseg
    %On extrait un segment de taille M du signal
    seg(i,:)=input(1+(i-1)*n:(i-1)*n+M);
    %On cherche le polynome de degré K qui correspond
    P=polyfit(index_y,seg(i,:),K);
    y(i,:)=polyval(P,index_y);
 
end
 
%On met tous les segments sous forme d'un vecteur
B=reshape(y',1,size(y,1).*size(y,2));
 
A=[];
for i=1:Nseg-1
%D'apres l'article, on calcule la valeur de Yc sur les 
%abscisses correspondant aux intersections de chaque segment
 
Yc(i,index_yc)=y(i,index_yc+n).*(1-(index_yc-1)/n)+y(i+1,index_yc).*(index_yc-1)/n;
end
 
%Mise des Yc sous forme d'un seul vecteur
A=reshape(Yc',1,size(Yc,1).*size(Yc,2))';

En relisant mon code j'ai changé 2/3 ça ne compile plus. J'ai des problèmes de taille. Merci de votre aide, mais je pense qu'en fait le soucis vient que je sais pas bien compter....
Il faut bien calculer la valeur de Nseg: le nombre de segment de taille 2n+1 que l'on peut mettre dans un axe de longueur length(input) et dont chaque segment recouvrement la moitié (n+1 points) de son voisin de droite et est recouvert par la moitié de son voisin de gauche (n+1 points).

**tachmou** · 16/08/2013, 10h08

Ah, c'est déjà bien plus clair !

Bon, je me suis permis d'apporter quelques modifications

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
clc;
close all;
clear all;

% Input
input = [];
load 'input.mat'
% input est de taille 800x1

n=4;

% Taille d'un segment
M=2*n+1;

% Nombre de segments de taille M qui recouvrent le signal
% et qui ont n+1 points en commun
Nseg=floor((length(input)-(n+1))/n)-1;

% Degre du polynome y
K=2;

% Index de l'axe des abscisse de chaque polynôme y
index_y=1:M;

% Index de l'axe des abcsisse de chaque Yc
index_yc=1:n+1;


B=[];

for i=1:Nseg
    % On extrait un segment de taille M du signal
    seg(i,:)=input(1+(i-1)*n:(i-1)*n+M);

    % On cherche le polynome de degré K qui correspond
    P=polyfit(index_y,seg(i,:),K);
    y(i,:)=polyval(P,index_y);
end

%On met tous les segments sous forme d'un vecteur
B=reshape(y',1,size(y,1).*size(y,2));

A=[];

for i=1:Nseg-1
    % D'apres l'article, on calcule la valeur de Yc sur les 
    % abscisses correspondant aux intersections de chaque segment
    % y(c)(l)=(1-(l-1)/n)*y(i)(l+n) + (l-1)/n*y(i+1)(l), l=1,2,..,n+1

    Yc(i,index_yc)=y(i,index_yc+n).*(1-(index_yc-1)/n)+y(i+1,index_yc).*(index_yc-1)/n;
end

%Mise des Yc sous forme d'un seul vecteur
A=reshape(Yc',1,size(Yc,1).*size(Yc,2))';

(==> pas testé, je n'ai pas matlab sous la main en ce moment)

Pas la peine de mettre sous forme de vecteur : ta première méthode était bonne... quand je te disais de travailler sur des séries à une dimension, je voulais parler de ton "input" car à la base j'avais peur que justement tu travaillais sur plusieurs séries en même temps ^^'

Quelles sont les erreurs dont tu parles ?

Par ailleurs, maintenant je peux comprendre ta remarque du début :

Ce qui m’étonne (voir m’énerve) c'est qu'au final le polynôme Yc que j'obtiens est plus long que le signal de base... Alors que normalement, avec la division euclidienne on devrait tomber sur une taille au moins plus petite.

==> Je ne sais pas d'où tu sort ce "normalement", mais pour moi cette vérité est ... fausse ! Si on fait un calcul "à la louche" : Nseg*(n+1) > Nseg*n ! (Attention, ceci n'est vrai qu'à condition d'avoir beaucoup de découpe, bien sûr... sinon ça peut être effectivement plus petit dans certains cas bien particuliers : quand le reste de la division euclidienne est plus grande que Nseg !) En effet, dans le signal de départ, tes signaux se recouvrent. A la fin, ton signal ne se recouvre pas mais se juxtapose... les "1" se "dédoublent" en quelque sorte...

Du coup, ton code de départ était peut-être juste

... Content d'avoir passé quelques jours pour rien à cause de moi ?

N'empêche, comme quoi, ça vaut le coup de prendre le temps de rendre son code et ses explications le plus didactique possible, surtout quand on demande de l'aide, même si on a l'impression qu'on prend les gens pour des c*** si on fait ça : c'est tout le contraire en fait, donc il ne faut pas hésiter ^^ ! D'ailleurs, bien des fois, en détaillant tranquillement tout, on trouve l'erreur soi-même très rapidement du coup... et si on ne la trouve pas, elle devient évidente pour un individu lambda ==> tout bénéf!

J'espère que ça répond à tes questions

ps: je ne suis pas très fan des ".*", même si ça allège énormément le code, je trouve que ça le rend plus difficilement lisible... du moins c'est plus difficile de trouver des erreurs de calculs avec... mais là, ça me semble juste.

**darkwhite** · 19/08/2013, 18h59

Merci infiniment pour ton aide. Je rentre de voyage donc j'étudierais ta solution en profondeur un peu plus tard.

==> Je ne sais pas d'où tu sort ce "normalement", mais pour moi cette vérité est ... fausse ! Si on fait un calcul "à la louche" : Nseg*(n+1) > Nseg*n ! (Attention, ceci n'est vrai qu'à condition d'avoir beaucoup de découpe, bien sûr... sinon ça peut être effectivement plus petit dans certains cas bien particuliers : quand le reste de la division euclidienne est plus grande que Nseg !) En effet, dans le signal de départ, tes signaux se recouvrent. A la fin, ton signal ne se recouvre pas mais se juxtapose... les "1" se "dédoublent" en quelque sorte...

Je crois que mon soucis est là, je voulais me débarrasser des points "en trop" ou "compté deux fois". Car quand je disais que je ne comprenais pas pourquoi le signal des intersections était plus grand que le signal de départ c'était par vision géométrique sur le segment...
Au fond, qu'il y ait des points que l'on compte deux fois ne me dérange pas plus que ça, le problème c'est que je dois soustraire le polynôme obtenu au signal de départ. Donc s'il ne sont pas superposé, il n'y a plus beaucoup d’intérêt à faire un recouvrement fin comme celui là ...
Je me trompe encore sur les coupures à operer dans le signal pour avoir une superposition (par rapport à l'axe des abscisses bien sur).

Merci encore

**tachmou** · 20/08/2013, 09h24

Bonjour,

Envoyé par darkwhite

Je crois que mon soucis est là, je voulais me débarrasser des points "en trop" ou "compté deux fois". Car quand je disais que je ne comprenais pas pourquoi le signal des intersections était plus grand que le signal de départ c'était par vision géométrique sur le segment...

Je ne comprends pas trop ce passage hmm.. C'est bien par vision géométrique, il n'y a pas de soucis... Les Nseg signaux que tu obtiens à la fin sont de taille (n+1) alors que ceux utilisés au début sont de taille (n) en comparaison, à quelques détails près... donc, selon les cas, tu peux avoir soit plus grand, soit plus petit. Dans le cas que tu as cité en dernier, le signal obtenu à la fin est plus grand que celui de départ : imaginons que Nseg = 300, et n = 3 alors tailleSignalFinMisBoutABout = Nseg*(n+1) = 1200 > 900 = Nseg*n = tailleSignalDépart

Envoyé par darkwhite

Au fond, qu'il y ait des points que l'on compte deux fois ne me dérange pas plus que ça, le problème c'est que je dois soustraire le polynôme obtenu au signal de départ.
Donc s'il ne sont pas superposé, il n'y a plus beaucoup d’intérêt à faire un recouvrement fin comme celui là ...
Je me trompe encore sur les coupures à operer dans le signal pour avoir une superposition (par rapport à l'axe des abscisses bien sur).

Alors là, je ne vois pas du tout de quoi tu parles... en tout cas je n'arrive pas à faire le lien avec ton premier post et le bout d'article que tu as cité... bon, après, je ne suis pas un boss du signal, donc il y a peut-être une chose "évidente" pour un pro du signal que je ne sais pas...

Cordialement,

**darkwhite** · 20/08/2013, 23h13

Envoyé par tachmou

Bonjour,

Je ne comprends pas trop ce passage hmm.. C'est bien par vision géométrique, il n'y a pas de soucis... Les Nseg signaux que tu obtiens à la fin sont de taille (n+1) alors que ceux utilisés au début sont de taille (n) en comparaison, à quelques détails près... donc, selon les cas, tu peux avoir soit plus grand, soit plus petit. Dans le cas que tu as cité en dernier, le signal obtenu à la fin est plus grand que celui de départ : imaginons que Nseg = 300, et n = 3 alors tailleSignalFinMisBoutABout = Nseg*(n+1) = 1200 > 900 = Nseg*n = tailleSignalDépart

Alors là, je ne vois pas du tout de quoi tu parles... en tout cas je n'arrive pas à faire le lien avec ton premier post et le bout d'article que tu as cité... bon, après, je ne suis pas un boss du signal, donc il y a peut-être une chose "évidente" pour un pro du signal que je ne sais pas...

Cordialement,

Je dois me tromper quelque part mais la vision géométrique dont je parlais c'était que l'on calcul une quantité sur les intersections des segments (sur mon image moche, les zones hachuré en gris. Et donc sur la reunion des zones grises on pourrait relier point a point le signal initial et le fitting polynomial (ou le soustraire).

**tachmou** · 22/08/2013, 10h14

Bon, on se perd là j'ai l'impression.

Ton schéma manque de clarté je trouve, je te conseille d'en faire un plutôt en cascade. Imaginons que je symbolise un bout de signal de taille n par "===" et un bout du signal de taille 1 par "-" :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
-------------------
===-===
   ===-===
      ===-===
         ===-===
            ===-===

Sur mon exemple, mon "input" est de taille 19.
n = 3
M = 2*n+1 = 7
Nseg = floor((19-(n+1))/n) = floor((19-4)/3) = 5 (qui aurait donné le même résultat avec 20 et 21, mais avec 22 on aurait obtenu 6, ce qui est normal)
A la fin, tes bouts de signaux Yc sont de taille n+1 = 4. Il y en a en tout Nseg = 5. Donc si tu les met bout à bout, tu obtiens un signal complet de 20 élément > 19. CQFD

Tu observeras que si ton signal de départ vaut 19, 20 ou 21, tu obtiens quoi qu'il en soit des bouts de signaux Yc qui, si tu les colles bout à bout, seront de taille 20. Donc selon la taille de ton signal de départ, tu peux obtenir un signal soit plus grand, soit plus petit. Tout dépend de la taille de ton input, de n et de Nseg.

Après, "coller" les signaux bout à bout n'est peut-être pas la bonne méthode (certainement même). Mais je ne connais pas l'article et je ne suis pas un pro du signal : donc à toi de voir quelle est la "bonne" méthode.

Cordialement,

ps: à quoi est censé servir cet algorithme ? J'ai l'impression que c'est beaucoup de calcul pour pas grand chose. Un simple filtrage avec un schéma pondéré centré pourrait suffire si tu cherches juste à lisser le signal...

**darkwhite** · 22/08/2013, 17h02

Merci infiniment, je pense que mes petits synapses commencent à admettre mon erreur.
Effectivement c'est beaucoup de prise de tête pour pas grand chose. Le but etait en fait d'enlever ce que cet article anglais appelle "a sinusoidal trend" à une série qui représente le nombre de taches solaire. Ils comparait ça à la méthode bourrine de faire un passe bas brut en fourrier en passant par une fft tronquée.
En fait le problème est plus devenu une question casse tête qu'indispensable...
Si vous êtes toujours interessé et par curiosité, l'article complet :
http://iopscience.iop.org/1742-5468/..._02_P02066.pdf

Detrending algorithme et nombre de points

Signal

Discussions similaires

Partager

Partager