Deconvolução de sinais 1D embaçados por um núcleo gaussiano

Convolvi um sinal aleatório com um Gaussiano e adicionei ruído (nesse caso, ruído de Poisson) para gerar um sinal barulhento. Agora eu gostaria de desconvolver esse sinal barulhento para extrair o sinal original usando o mesmo gaussiano.

O problema é que preciso de um código que faça o trabalho de deconvolução em 1D. (Eu já encontrei alguns em 2D, mas meu objetivo principal é 1D).

Você pode me sugerir alguns pacotes ou programas capazes de fazer isso? (De preferência em MATLAB)

Obrigado antecipadamente pela ajuda.

Eu expliquei uma vez no StackOverflow .

---

Seu sinal pode ser representado como um vetor e a convolução é a multiplicação com uma matriz diagonal N (onde N é o comprimento do filtro). Estou assumindo, por uma questão de resposta, que o filtro é muito menor que o sinal

Por exemplo:

Seu vetor / sinal é:

V1
V2
...
Vn

Seu filtro (elemento de convolução) é:

  [b1 b2 b3];

Portanto, a matriz é nxn: (Seja chamada A):

[b2 b3 0  0  0  0.... 0]
[b1 b2 b3 0  0  0.... 0]
[0  b1 b2 b3 0  0.... 0]
.....
[0  0  0  0  0  0...b2 b3]

Convolução é:

A*v;

E a desconvolução é

A^(-1) * ( A) * v;

Obviamente, em alguns casos a desconvolução não é possível. Estes são os casos em que você possui A. singular. Mesmo matrizes que não são singulares, mas próximas de serem singulares, podem ser problemáticas, pois terão um grande erro numérico. Você pode estimar calculando o número da condição da matriz.

Se A tiver uma condição baixa, você poderá calcular o inverso e aplicá-lo no resultado.

---

Agora, vamos ver alguns exemplos no Matlab:

Primeiro, criei uma função que calcula a matriz de convolução.

function A = GetConvolutionMatrix(b,numA)
    A = zeros(numA,numA);
    vec = [b  zeros(1,numA-numel(b))];
    for i=1:size(A,1)
        A(i,:) = circshift(vec,[1 i]);
    end
end

Agora, vamos tentar ver o que acontece com diferentes kernels:

    b = [1 1 1];
    A = GetConvolutionMatrix(b,10);
    disp(cond(A));

O número da condição é:

 7.8541

Este é problemático, como esperado. Após a média, é difícil recuperar o sinal original.

Agora, vamos tentar uma média mais branda:

b = [0.1 0.8 0.1];
A = GetConvolutionMatrix(b,10);
disp(cond(A));

O resultado é:

1.6667

Isso vai bem com a nossa intuição. É muito mais fácil reverter a média leve do sinal original.

Também podemos ver como a matriz inversa se parece:

 figure;imagesc(inv(A));

Aqui está uma linha da matriz:

  0.0003   -0.0026    0.0208   -0.1640    1.2910   -0.1640    0.0208   -0.0026    0.0003   -0.0001

Podemos ver que a maior parte da energia em cada linha está concentrada em 3-5 coeficientes ao redor do centro. Portanto, para deconvolver, podemos simplesmente convencer o sinal novamente com esta aproximação:

   [0.0208   -0.1640    1.2910   -0.1640    0.0208]

Este kernel parece interessante! É um operador de afiação. Nossa intuição está correta, a nitidez cancela a desfocagem.

Se você adicionou ruído aleatório, não pode obter o sinal original ... Você pode tentar separar os sinais no domínio da frequência (se o ruído e o sinal forem de frequências diferentes). Mas parece que o que você está procurando é um filtro Wiener .

Eu acho que isso ainda é um problema em aberto.

Existem inúmeros trabalhos de pesquisa que tentam recuperar o sinal original da melhor maneira possível.

Uma abordagem clássica é através dos métodos baseados no Wavelet .

Existem também abordagens de dicionário como esta .

Você pode ter uma visão mais aprofundada do problema seguindo a pesquisa de David L. Donho, Michael Elad, Alfred M. Bruckstein etc.

Se entendi o problema corretamente, podemos formalizá-lo da seguinte maneira:

Nós temos um modelo de sinal,

$y = Hx + \eta$

$y$$H$$\eta$$x$

$\eta$

Não trabalhei na recuperação de sinal sob o ruído de Poisson, mas pesquisei e achei este documento útil. Abordagens semelhantes nesse contexto podem ser úteis para esse problema.

Sabe-se que a desconvolução de dados ruidosos é um problema incorreto, uma vez que o ruído é arbitrariamente ampliado no sinal reconstruído. Portanto, é necessário um método de regularização para estabilizar a solução. Aqui, você pode encontrar um pacote MATLAB que resolve esse problema implementando o algoritmo de regularização de Tikhonov:

https://github.com/soheil-soltani/TranKin .

Eu irei para o começo da questão. Existem funções de desconvolução no MATLAB que são usadas para aplicativos de processamento de imagem. No entanto, você também pode usar essas funções para sinais 1D. Por exemplo,

% a random signal
sig_clean = zeros(1,200); 
sig_clean(80:100)=100;

figure
subplot(1,3,1)
plot(sig_clean,'b-.','LineWidth',2)
legend('Clean Signal')

% convolve it with a gaussian
x=1:30;
h = exp(-(x-15).^2/20); h=h/sum(h);
sig_noisy = conv(sig_clean,h,'same');

% and add noise
sig_noisy = awgn(sig_noisy,0,'measured');

subplot(1,3,2)
plot(sig_noisy,'r')
hold on, plot(sig_clean,'b-.','LineWidth',3)
legend('Blurred and noise added signal','Clean Signal')

( sig_noisy = sig_clean * h + noise) Então porque não desconvolver o sinal de saída com a hfunção e obter o (quase) sinal de entrada. Estou usando a deconvolução Wiener aqui

sig_deconvolved=deconvwnr(sig_noisy,h,1);

subplot(1,3,3)
plot(sig_noisy,'r')
hold on, plot(sig_clean,'b-.','LineWidth',2)
hold on, plot(sig_deconvolved,'k--','LineWidth',2)
legend('Blurred and noise added signal','Clean Signal','Deconvolved Signal')

Como alternativa, se você não conhece a hfunção, mas conhece a entrada e a saída, desta vez, por que não desconvolver o sinal de entrada com a saída que fornecerá a h^-1função. Então você pode usá-lo como um filtro para filtrar o sinal barulhento. ( sig_clean = sig_noisy * h^-1)

h_inv=deconvwnr(sig_clean,sig_noisy,1);

figure;
subplot(1,2,1)
plot(h_inv)
legend('h^-^1')


sig_filtered=conv(sig_noisy,h_inv,'same');
subplot(1,2,2)
plot(sig_noisy,'r')
hold on, plot(sig_clean,'b-.','LineWidth',2)
hold on, plot(sig_filtered,'k--','LineWidth',2)
legend('Blurred and noise added signal','Clean Signal','Filtered Signal')

Espero que ajude.

Convolução é a multiplicação e soma de dois sinais, um sobre o outro. Estou falando de dois sinais determinísticos. Se você deseja desconvolver uma da outra, isso corresponde à solução do sistema de equações. Como você deve saber, o sistema de equações nem sempre é solucionável. O sistema de equações pode ser superdeterminado, subdeterminado ou exatamente solucionável.

Caso você adicione algum ruído, você perde algumas informações e não pode recuperá-las. O que você pode fazer é novamente resolver o sistema linear de equações, considerando o fato de que cada coeficiente é adicionado a um termo de ruído. Ou, como você pode ver em outra resposta à sua pergunta, você pode primeiro estimar o sinal original do sinal barulhento e depois tentar resolver o sistema de equações.

É importante notar que o ruído é adicionado aos coeficientes multiplicados e somados. Portanto, pode até ser o caso de seu sistema de equações eventualmente não ser solucionável de maneira única. Para ter certeza de que é solucionável de forma exclusiva, sua matriz de coeficientes deve ser quadrada e de classificação completa.

Isso seria difícil de fazer. Convolução com um Gaussiano é equivalente à multiplicação com uma Transformada de Fourier do Gaussiano no domínio da frequência. Por acaso também é um gaussiano, em essência, é um filtro passa-baixo e é realmente eficaz nisso. Depois de adicionar ruído, todas as informações que estão na "faixa de parada" do Gaussian são destruídas. Não há como recuperar isso.

A desconvolução está multiplicando-se essencialmente com o inverso da resposta em frequência. Aqui está o problema: o inverso da resposta em frequência fica muito, muito grande onde o gaussiano original é muito pequeno. Nessas frequências, você basicamente amplificaria o ruído em grandes quantidades. Mesmo se tudo fosse completamente sem ruído, você provavelmente teria problemas numéricos.

Abordagens

Existem muitos métodos para a Deconvolução (ou seja, o operador de degradação é linear e invariante no tempo / espaço) por aí.
Todos eles tentam lidar com o fato de que o problema está mal resolvido em muitos casos.

Métodos melhores são aqueles que adicionam alguma regularização ao modelo dos dados a serem restaurados.
Pode ser modelos estatísticos (Priors) ou qualquer conhecimento.
Para imagens, um bom modelo é suave ou esparso dos gradientes.

Mas, para obter a resposta, será necessária uma abordagem paramétrica simples - Minimizar o erro dos mínimos quadrados entre os dados restaurados no modelo e as medições.

Modelo

O modelo de mínimos quadrados é simples.
A função objetivo em função dos dados é dada por:

$$f \left( x \right) = \frac{1}{2} \left\| h \ast x - y \right\|_{2}^{2}$$

O problema de otimização é dado por:

$$\arg \min_{x} f \left( x \right) = \arg \min_{x} \frac{1}{2} \left\| h \ast x - y \right\|_{2}^{2}$$

$x$$h$$y$
$x \in \mathbb{R}^{n}$$h \in \mathbb{R}^{k}$$y \in \mathbb{R}^{m}$$m = n - k + 1$

Esta é uma operação linear em espaço finito, portanto, pode ser escrita usando um formulário de matriz:

$$\arg \min_{x} f \left( x \right) = \arg \min_{x} \frac{1}{2} \left\| H x - y \right\|_{2}^{2}$$

Onde $H \in \mathbb{R}^{m \times n}$ é a matriz de convolução.

Solução

A solução dos mínimos quadrados é dada por:

$$\hat{x} = \left( {H}^{T} H \right)^{-1} {H}^{T} y$$

Como pode ser visto, requer uma inversão da matriz.
A capacidade de resolver isso adequadamente depende do número da condição do operador${H}^{T} H$ que obedece $\operatorname{cond} \left( H \right) = \sqrt{\operatorname{cond} \left( {H}^{T} H \right)}$.

Análise do número da condição

O que há por trás desse número de condição?
Pode-se responder usando Álgebra Linear.
Mas uma abordagem mais intuitiva, na minha opinião, seria pensar no domínio da frequência.

Basicamente, o operador de degradação atenua a energia, geralmente, de alta frequência.
Agora, como na frequência isso é basicamente uma multiplicação por elementos, alguém poderia dizer que a maneira mais fácil de invertê-la é a divisão por elemento pelo filtro inverso.
Bem, é o que é feito acima.
O problema surge com os casos em que o filtro atenua a energia praticamente em zero. Então temos problemas reais ...
Isso é basicamente o que o Número da Condição nos diz, quão difícil algumas frequências foram atenuadas em relação a outras.

Acima, era possível ver o número da condição (usando unidades [dB]) como uma função do parâmetro STD do filtro gaussiano.
Como esperado, quanto maior a DST, pior o número da condição, pois maior a DST significa LPF mais forte (os valores que caem no final são questões numéricas).

Solução Numérica

O conjunto de Gaussian Blur Kernel foi criado.

Os parâmetros são $n = 300$, $k = 31$ e $m = 270$.
Os dados são aleatórios e nenhum ruído foi adicionado.

No MATLAB, o Sistema Linear foi resolvido usando o pinv()Pseudo Inverso baseado em SVD e o \operador.

Como se pode ver, usando o SVD, a solução é muito menos sensível conforme o esperado.

Por que existe um erro?
Procurando uma solução (para a maior DST):

Como se pode ver, o sinal é restaurado muito bem, exceto no início e no fim.
Isso se deve ao uso da Valid Convolution, que nos diz pouco sobre essas amostras.

Ruído

Se adicionássemos ruído, as coisas pareceriam diferentes!
A razão pela qual os resultados eram bons antes se deve ao fato do MATLAB poder manipular o DR dos dados e resolver as equações, mesmo que eles tenham um número de condição grande.

Porém, um grande número de condições significa que o filtro inverso amplifica fortemente (para reverter a atenuação forte) algumas frequências.
Quando estes contêm ruído, significa que o ruído será amplificado e a restauração será ruim.

Como se pode ver acima, agora a reconstrução não funcionará.

Sumário

Se alguém souber exatamente o Operador de Degradação e o SNR for muito bom, métodos simples de desconvolução funcionarão.
A principal questão da desconvolução é a intensidade com que o operador de degradação atenua as frequências.
Quanto mais atenua, mais SNR é necessário para restaurar (esta é basicamente a ideia por trás do Wiener Filter ).
As frequências que foram definidas como zero não podem ser restauradas!

Na prática, a fim de obter resultados estáveis, deve adicionar alguns antecedentes.

O código está disponível no meu repositório StackExchange Signal Processing Q2969 GitHub .

Em geral, um método para lidar com a questão que generaliza substancialmente um problema de extração de dois ou mais componentes é pegar os espectros G¹, G² ⋯, Gⁿ dos sinais # 1, # 2, ..., #n, tabular o total quadrado Γ (ν) = | G¹ (ν) | ² + | G² (ν) | ² + ⋯ + | Gⁿ (ν) | ² em cada frequência ν e normalize G₁ (ν) ≡ G¹ (ν) * / Γ (ν), G₂ (ν) ≡ G² (ν) * / Γ (ν), ..., G_n (ν) ≡ Gⁿ (ν) * / Γ (ν). O problema com a falta de definição e o ruído corresponde ao fato de que ν (ν) ~ 0 é possível para algumas frequências ν. Para lidar com isso, adicione outro "sinal" para extrair G⁰ (ν) = constante - o sinal de "ruído". Agora Γ (ν) será estritamente limitado abaixo. Isso está quase certamente relacionado à regularização de Tikhonov, mas nunca encontrei ou estabeleci nenhum resultado de equivalência ou outra correspondência. É mais simples, mais direto e intuitivo.

Como alternativa, você pode tratar os Gs como vetores equipados com um produto interno adequado, por exemplo, «G, G '» ∫ ∫ G (ν) * G' (ν) dν e tomar (G₀, G₁, ⋯, G_n) como um duplo (por exemplo, o inverso generalizado) de (G⁰, G¹, ⋯, Gⁿ) - assumindo, é claro, que os vetores componentes são linearmente independentes.

Para a desconvolução gaussiana, seria estabelecido n = 1, G⁰ = sinal de "ruído" e G¹ = sinal de "Gaussian".

A resposta fornecida por Andrey Rubshtein falhará miseravelmente na presença de ruído, pois o problema descrito é muito sensível a ruído e erros de modelagem. É uma boa idéia construir uma matriz de convolução, mas o uso da regularização na inversão é uma necessidade absoluta em um problema como esse. Um método de regularização muito simples e direto (embora computacionalmente caro) é a Decomposição de Valor Singular Truncado que (TSVD). Métodos como Regularização de Tikhonov e vale a pena conferir. A regularização de Tikhonov (e sua forma geral) possui uma forma empilhada muito elegante e fácil de implementar no Matlab. Confira o livro: Problemas inversos lineares e não lineares com aplicações práticas de Samuli Siltanen e Jennifer Mueller. Total de Variação

Na verdade, a questão não está clara. Mas as respostas carificaram o que você pediu. Você pode construir um sistema de equações algébricas lineares, como algumas pessoas aconselham, isso é correto, mas a matriz construída no sinal conhecido é chamada de mal condicionada. Isso significa que quando você tenta invertê-lo, os erros de truncamento eliminam a solução e você recebe números aleatórios no resultado. A abordagem comum é restrita ao extremo. Você minimiza a norma da solução || x || com restrição || Ax - y || <delta. Então você está procurando um x com a menor norma que não permita que a diferença entre Ax e y seja grande. É muito simples que você precise adicionar o chamado parâmetro de regularização na principal diagonal da matriz obtida na aplicação dos mínimos quadrados. É chamado regularização de Tikhonov. Eu tenho amostras de codificação que fazem isso,