Valid Shape FFT Deconvolution

11 visualizaciones (últimos 30 días)
Benjamin Hezrony
Benjamin Hezrony el 6 de Oct. de 2023
Comentada: Benjamin Hezrony el 18 de Oct. de 2023
Hi MathWorks community,
I am trying to wrap my head around FFT convolution and deconvolution, and the effect of different "shape" (full, valid, and same). I have sample code displaying my problem. I'd like to reconstruct the input, x, in both deconvolutions. I got it for full conv/deconv, but not for valid conv/deconv. I didn't even try it for same conv/deconv yet since I don't even know if it is possible for valid conv/deconv. I hope it is possible. I may be missing some theory!
Thank you for any help,
~BH

Respuesta aceptada

Paul
Paul el 7 de Oct. de 2023
Hi Benjamin,
Here's an approach for the 'valid' convolution.
Problem parameters:
x = 1:7; nx = numel(x);
h = 1:3; nh = numel(h);
nyf = nx + nh - 1;
nyv = nx - nh + 1;
Define the window function such that conv(x,h,'valid) = w*conv(x,h,'full') with appropriate zero padding as will be shown
d = nh - 1;
w = [zeros(1,d) ones(1,nyv) zeros(1,nyf - nyv - d)];
nw = numel(w);
Verify
yf = conv(x,h)
yf = 1×9
1 4 10 16 22 28 34 32 21
w.*yf
ans = 1×9
0 0 10 16 22 28 34 0 0
conv(x,h,'valid')
ans = 1×5
10 16 22 28 34
Now, just like we could compute yf as
yf = ifft(fft(x,nw).*fft(h,nw)) % same results as above
yf = 1×9
1.0000 4.0000 10.0000 16.0000 22.0000 28.0000 34.0000 32.0000 21.0000
We can compute yv (the zero-padded form of the valid convolution) as
yv = ifft(cconv(fft(w),fft(x,nyf).*fft(h,nyf),nw)/nw) % same as above (*)
yv =
0.0000 + 0.0000i 0.0000 + 0.0000i 10.0000 + 0.0000i 16.0000 + 0.0000i 22.0000 - 0.0000i 28.0000 + 0.0000i 34.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 - 0.0000i
Verify the imaginary part is rounding error and get rid of it
max(abs(imag(yv)))
ans = 2.3614e-15
yv = real(yv)
yv = 1×9
0 0.0000 10.0000 16.0000 22.0000 28.0000 34.0000 -0.0000 0.0000
From which conv(x,h,'valid') can be extracted with appropriate indexing.
If we want to recover x starting from fft(yv), it may be possible to invert (*) to get to x. But I couldn't see a way to do that. It seems clear that
fft(yv) = fft(w.*conv(x,h)) = fft(w.*ifft(fft(x,nw).*fft(h,nw)))
but I don't see a way to isolate an expression for fft(x), which could then be ifft'd.
Will be interested to see if anyone comes up with a solution.

Más respuestas (1)

Bruno Luong
Bruno Luong el 7 de Oct. de 2023
Editada: Bruno Luong el 7 de Oct. de 2023
You just do something that is obviously impossible. By FFT with truncation
X = fft(x,L)
You don't use the n-1 last elements so The (n-1) last elements of X is not used to get y.
There is then no way to get back X from Y and h2/H2
%% Valid conv and deconv:
x = 1:7;
h = 1:3;
m = length(x);
n = length(h);
L = m-n+1
L = 5
X = fft(x,L);
H = fft(h,L);
Y = X.*H;
y = ifft(Y,L)
y = 1×5
23 19 10 16 22
% This x3/h will give the same vector in Fourier domain as x/h
x3 = x; x3([6 7]) = 0;
X3 = fft(x3,L);
Y3 = X3.*H;
y3 = ifft(Y3,L)
y3 = 1×5
23 19 10 16 22
AFAIK the discrete FFT theorem applies for circular convolution. Only FULL convolution is can be convert to FT by Padding 0s and make pb circular without interfering the head and the tail of the original data.
One cannot do that with 'SAME' and 'VALID' convolution. This is impossible, period.
Another way to conceive that only FULL convolution is invertible (interm of x as unknown) is to write down the convolution matrix (in term of flip(h) as band diagonal). You'll see that the simple analysis is matrix size that only FULL convolution makes sense of possibme inversion.
But SAME and VALID convolution are simply truncations of the FULL convolution (that is how @Paul did in his reply).
  1 comentario
Benjamin Hezrony
Benjamin Hezrony el 18 de Oct. de 2023
@Paul, @Bruno Luong I wish I could formaly accept both of your answers! With your input and some further digging, I've come to understand why what I had been trying was impossible. Thank you both!
~BH

Iniciar sesión para comentar.

Categorías

Más información sobre Fourier Analysis and Filtering en Help Center y File Exchange.

Productos


Versión

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by