J'ai convolué un signal aléatoire avec un bruit gaussien et ajouté (bruit de Poisson dans ce cas) pour générer un signal bruyant. Maintenant, je voudrais déconvoluer ce signal bruyant pour extraire le signal d'origine en utilisant le même gaussien.
Le problème est que j'ai besoin d'un code qui fait le travail de déconvolution en 1D. (J'en ai déjà trouvé en 2D mais mon objectif principal est le 1D).
Pouvez-vous me suggérer des packages ou des programmes capables de le faire? (De préférence dans MATLAB)
Merci d'avance pour l'aide.
Réponses:
Je l'ai expliqué une fois sur StackOverflow .
Votre signal peut être représenté comme un vecteur et la convolution est une multiplication avec une matrice N-diagonale (où N est la longueur du filtre). Je suppose pour la réponse que le filtre est beaucoup plus petit que le signal
Par exemple:
Votre vecteur / signal est:
Votre filtre (élément convolving) est:
Donc, la matrice est nxn: (Qu'on l'appelle A):
La convolution est:
Et la déconvolution est
De toute évidence, dans certains cas, la déconvolution n'est pas possible. Ce sont les cas où vous avez un A. singulier. Même les matrices qui ne sont pas singulières, mais proches de l'être, peuvent être problématiques, car elles auront une grande erreur numérique. Vous pouvez l'estimer en calculant le numéro de condition de la matrice.
Si A est faible, vous pouvez calculer l'inverse et l'appliquer sur le résultat.
Voyons maintenant quelques exemples dans Matlab:
Tout d'abord, j'ai créé une fonction qui calcule la matrice de convolution.
Maintenant, essayons de voir ce qui se passe avec différents noyaux:
Le numéro de condition est:
Celui-ci est problématique, comme prévu. Après la moyenne, il est difficile de récupérer le signal d'origine.
Maintenant, essayons une moyenne plus douce:
Le résultat est:
Cela va bien avec notre intuition, une moyenne douce du signal d'origine est beaucoup plus facile à inverser.
Nous pouvons également voir à quoi ressemble la matrice inverse:
Voici une ligne de la matrice:
Nous pouvons voir que la plupart de l'énergie dans chaque ligne est concentrée en 3-5 coefficients autour du centre. Par conséquent, afin de déconvoluer, nous pouvons simplement convoluer à nouveau le signal avec cette approximation:
Ce noyau semble intéressant! C'est un opérateur d'affûtage. Notre intuition est correcte, la netteté annule le flou.
la source
Si vous avez ajouté du bruit aléatoire, vous ne pouvez pas obtenir le signal d'origine ... Vous pouvez essayer de séparer les signaux dans le domaine fréquentiel (si le bruit et le signal sont de fréquences différentes). Mais il semble que ce que vous recherchez soit un filtre de Wiener .
la source
Je pense que c'est toujours un problème ouvert.
Il existe de nombreux articles de recherche qui tentent de récupérer le signal d'origine du mieux qu'ils peuvent.
Une approche classique consiste à utiliser des méthodes basées sur les ondelettes .
Il existe également des approches de dictionnaire comme celle- ci.
Vous pouvez obtenir une vue plus approfondie du problème en suivant les recherches effectuées par David L. Donho, Michael Elad, Alfred M. Bruckstein, etc.
la source
Si j'ai bien compris le problème, nous pouvons formaliser le problème comme suit:
Nous avons un modèle de signal,
Je n'ai pas travaillé sur la récupération du signal sous le bruit de Poisson, mais j'ai googlé et j'ai trouvé que cet article pouvait être utile. Des approches similaires dans ce contexte peuvent être utiles pour ce problème.
la source
La déconvolution d'une donnée bruyante est connue pour être un problème mal posé, car le bruit est amplifié arbitrairement dans le signal reconstruit. Par conséquent, une méthode de régularisation est nécessaire pour stabiliser la solution. Ici, vous pouvez trouver un package MATLAB qui résout ce problème en implémentant l'algorithme de régularisation de Tikhonov:
https://github.com/soheil-soltani/TranKin .
la source
J'irai au tout début de la question. MATLAB contient des fonctions de déconvolution qui sont utilisées pour les applications de traitement d'image. Cependant, vous pouvez également utiliser ces fonctions pour les signaux 1D. Par exemple,
(
sig_noisy = sig_clean * h + noise
) Alors pourquoi ne pas déconvoluer le signal de sortie avec lah
fonction et obtenir (presque) le signal d'entrée. J'utilise la déconvolution de Wiener iciAlternativement, si vous ne connaissez pas la
h
fonction, mais connaissez l'entrée et la sortie, cette fois pourquoi ne pas déconvoluer le signal d'entrée avec la sortie qui donnera lah^-1
fonction. Ensuite, vous pouvez l'utiliser comme filtre pour filtrer le signal bruyant. (sig_clean = sig_noisy * h^-1
)J'espère que ça aide.
la source
La convolution est la multiplication et la sommation de deux signaux l'un sur l'autre. Je parle de deux signaux déterministes. Si vous voulez déconvoluer l'un de l'autre, cela correspond à la solution du système d'équations. Comme vous le savez peut-être, les systèmes d'équations ne sont pas toujours résolubles. Le système d'équations peut être surdéterminé, sous-déterminé ou exactement résoluble.
Si vous ajoutez du bruit, vous perdez des informations et vous ne pouvez pas récupérer ces informations. Ce que vous pouvez faire, c'est à nouveau de résoudre le système linéaire d'équations en considérant le fait que chaque coefficient est ajouté un terme de bruit. Ou, comme vous pouvez le voir dans une autre réponse à votre question, vous voudrez peut-être d'abord estimer le signal d'origine à partir du signal bruyant, puis essayer de résoudre le système d'équations.
Il est important de noter que le bruit s'ajoute aux coefficients multipliés et additionnés. Par conséquent, il se pourrait même que votre système d'équation ne soit finalement pas uniquement résoluble. Pour être sûr qu'il est uniquement résoluble, votre matrice de coefficients doit être carrée et de rang complet.
la source
Ce serait difficile à faire. La convolution avec un gaussien équivaut à une multiplication avec une transformée de Fourier du gaussien dans le domaine fréquentiel. Il se trouve que c'est aussi un gaussien, c'est essentiellement un filtre passe-bas et vraiment efficace. Une fois que vous ajoutez du bruit, toutes les informations contenues dans la "bande d'arrêt" du gaussien sont détruites. Il n'y a aucun moyen de récupérer cela.
La déconvolution se multiplie essentiellement avec l'inverse de la réponse en fréquence. Voici le problème: l'inverse de la réponse en fréquence devient vraiment très grand lorsque le gaussien d'origine est très petit. À ces fréquences, vous amplifieriez essentiellement le bruit par des quantités énormes. Même si tout était complètement silencieux, vous rencontreriez très probablement des problèmes numériques.
la source
Approches
Il existe de nombreuses méthodes de déconvolution (à savoir, l'opérateur de dégradation est linéaire et invariant temps / espace).
Tous essaient de faire face au fait que le problème est empoisonné dans de nombreux cas.
Les meilleures méthodes sont celles qui ajoutent une certaine régularisation au modèle des données à restaurer.
Il peut s'agir de modèles statistiques (Prieurs) ou de toute connaissance.
Pour les images, un bon modèle est le lissage par morceaux ou la rareté des dégradés.
Mais pour la réponse, une approche paramétrique simple sera adoptée - Minimiser l'erreur des moindres carrés entre les données restaurées dans le modèle et les mesures.
Modèle
Le modèle des moindres carrés est simple.
La fonction objectif en fonction des données est donnée par:
Le problème d'optimisation est donné par:
Il s'agit d'une opération linéaire dans un espace fini et peut donc être écrite à l'aide d'une forme matricielle:
Solution
La solution des moindres carrés est donnée par:
Analyse du numéro de condition
Qu'est-ce qui se cache derrière ce numéro de condition?
On pourrait y répondre en utilisant l'algèbre linéaire.
Mais une approche plus intuitive, à mon avis, serait d'y penser dans le domaine fréquentiel.
Fondamentalement, l'opérateur de dégradation atténue l'énergie, généralement à haute fréquence.
Maintenant, comme en fréquence il s'agit fondamentalement d'une multiplication par élément, on dirait que la manière simple de l'inverser est la division par élément par le filtre inverse.
Eh bien, c'est ce qui est fait ci-dessus.
Le problème se pose avec les cas où le filtre atténue l'énergie pratiquement à zéro. Ensuite, nous avons de vrais problèmes ...
C'est essentiellement ce que le numéro de condition nous dit, à quel point certaines fréquences ont été atténuées par rapport à d'autres.
Au-dessus, on pouvait voir le numéro de condition (en utilisant des unités [dB]) en fonction du paramètre STD du filtre gaussien.
Comme prévu, plus le STD est élevé, plus le nombre de conditions est mauvais, car un STD plus élevé signifie un LPF plus fort (les valeurs qui diminuent à la fin sont des problèmes numériques).
Solution numérique
Ensemble de Gaussian Blur Kernel a été créé.
Dans MATLAB, le système linéaire a été résolu en utilisant
pinv()
Pseudo Inverse basé sur SVD et l'\
opérateur.Comme on peut le voir, en utilisant le SVD, la solution est beaucoup moins sensible que prévu.
Pourquoi y a-t-il une erreur?
Recherche d'une solution (pour la MST la plus élevée):
Comme on pouvait le voir, le signal est très bien restitué, sauf pour le début et la fin.
Cela est dû à l'utilisation de la convolution valide qui nous en dit peu sur ces échantillons.
Bruit
Si nous ajoutions du bruit, les choses auraient un aspect différent!
La raison pour laquelle les résultats étaient bons auparavant est due au fait que MATLAB pouvait gérer le DR des données et résoudre les équations même si leur nombre de conditions était élevé.
Mais un grand nombre de conditions signifie que le filtre inverse amplifie fortement (pour inverser la forte atténuation) certaines fréquences.
Lorsque ceux-ci contiennent du bruit, cela signifie que le bruit sera amplifié et que la restauration sera mauvaise.
Comme on pouvait le voir ci-dessus, maintenant la reconstruction ne fonctionnera pas.
Sommaire
Si l'on connaît exactement l'opérateur de dégradation et que le SNR est très bon, des méthodes de déconvolution simples fonctionneront.
Le problème principal de la déconvolution est de savoir à quel point l'opérateur de dégradation atténue les fréquences.
Plus il s'atténue, plus le SNR est nécessaire pour restaurer (c'est essentiellement l'idée derrière Wiener Filter ).
Les fréquences qui ont été mises à zéro ne peuvent pas être restaurées!
En pratique, pour avoir des résultats stables, il faut ajouter quelques priors.
Le code est disponible sur mon référentiel de traitement de signaux StackExchange Q2969 GitHub .
la source
En général, une méthode pour traiter le problème qui se généralise sensiblement à un problème d'extraction de deux composants ou plus consiste à prendre les spectres G¹, G² ⋯, Gⁿ des signaux # 1, # 2, ..., #n, tabuler le total carré Γ (ν) = | G¹ (ν) | ² + | G² (ν) | ² + ⋯ + | Gⁿ (ν) | ² à chaque fréquence ν, et normaliser G₁ (ν) ≡ G¹ (ν) * / Γ (ν), G₂ (ν) ≡ G² (ν) * / Γ (ν), ..., G_n (ν) ≡ Gⁿ (ν) * / Γ (ν). Le problème de mal-définition et de bruit correspond au fait que Γ (ν) ~ 0 est possible pour certaines fréquences ν. Pour gérer cela, ajoutez un autre "signal" pour extraire G⁰ (ν) = constante - le signal "bruit". Maintenant, Γ (ν) sera strictement limité ci-dessous. Ceci est presque certainement lié à la régularisation de Tikhonov, mais je n'ai jamais trouvé ni établi de résultat d'équivalence ou d'autre correspondance. C'est plus simple et plus direct et intuitif.
Alternativement, vous pouvez traiter les G comme des vecteurs équipés d'un produit intérieur approprié, par exemple «G, G '» ≡ ∫ G (ν) * G' (ν) dν, et prendre (G₀, G₁, ⋯, G_n) comme un dual (par exemple l'inverse généralisé) de (G⁰, G¹, ⋯, Gⁿ) - en supposant, bien sûr, que les vecteurs composants sont linéairement indépendants.
Pour la déconvolution gaussienne, on établirait n = 1, G⁰ = le signal "bruit" et G¹ = le signal "gaussien".
la source
La réponse fournie par Andrey Rubshtein échouera misérablement en présence de bruit, car le problème décrit est très sensible au bruit et aux erreurs de modélisation. C'est une bonne idée de construire une matrice de convolution, mais l'utilisation de la régularisation dans l'inversion est un must absolu dans un problème comme celui-ci. Une méthode de régularisation très simple et directe (bien que coûteuse en termes de calcul) est la décomposition de valeur singulière tronquée (TSVD). Des méthodes comme la régularisation de Tikhonov et la régularisation de la variation totalevalent le détour. La régularisation de Tikhonov (et sa forme générale) a une forme empilée très élégante qui est facile à implémenter dans Matlab. Consultez le livre: Problèmes inverses linéaires et non linéaires avec des applications pratiques par Samuli Siltanen et Jennifer Mueller.
la source
En fait, la question n'est pas claire. Mais les réponses ont confirmé ce que vous aviez demandé. Vous pouvez construire un système d'équations algébriques linéaires comme certains le conseillent, c'est exact, mais la matrice construite sur un signal connu est soi-disant mal conditionnée. Cela signifie que lorsque vous essayez de l'inverser, les erreurs de troncature tuent la solution et vous recevez des nombres aléatoires en résultat. L'approche courante est l'extrémum contraint. Vous minimisez la norme de solution || x || avec contrainte || Axe - y || <delta. Vous recherchez donc un tel x avec la plus petite norme qui ne permette pas à la différence entre Ax et y d'être grande. Il est très simple d'ajouter un soi-disant paramètre de régularisation sur la diagonale principale de la matrice obtenue en appliquant les moindres carrés. Cela s'appelle la régularisation de Tikhonov. J'ai des exemples de codage qui font ça,
la source