Reconciling convolution, fft and manual DFT results

Question

In Octave, I am playing with signal processing primitives, attempting to reproduce the convolution theorem in multiple ways: that convolution in the time domain is equivalent to point-wise multiplication in the frequency domain.

I consider three routes to reconstruct the original signal:

1. The fft and ifft functions
2. The conv function
3. A manually constructed DFT matrix.

I am attaching my working code, and the output, glad for inputs on where the bugs may be located.

N = 512;                  % number of points                                   
t = 0:N-1;                % [0,1,2,...,N-1]                                    
h = exp(-t);              % filter impulse reponse                             
H = fft(h);               % filter frequency response                          
x = (1+t) .* sin(sqrt(1+t)); % (input signal of our choice)                    
                                                                               
y1 = conv(x,h,"same");    % Direct convolution                                 
y2 = ifft(fft(x) .* H);   % FFT convolution                                    
                                                                               
T  = transpose(t) * t;                                                         
W  = exp(j * 2*pi/N * T); % DFT matrix                                         
y3 = (x * W .* H) * W/N;       % "Manual" convolution                          
                                                                               
lw = 2                                                                         
                                                                               
plot(t,                                                                        
    x,         ";orig;",                      "linewidth", lw+1,               
    y1,        ";conv;",   "linestyle", "--", "linewidth", lw,                 
    real(y2),  ";fft;",    "linestyle", ":",  "linewidth", lw,                 
    real(y3),  ";manual;", "linestyle", "-.", "linewidth", lw)                 
                                                                               
set(gca, "fontsize", 20, "linewidth", lw)                                      

1. In the first case (yellow), I am able to reconstruct the signal, but the scaling is not right (in fact the scaling is wrong in every case).
2. In the second case (red), it looks like the result is shifted, and half of the signal is lost.
3. In the third case (purple), I get something that's equivalent to fft but flipped horizontally.

Answer 1

Issues:

1. Your definition of the kernel is not right. conv expects the origin of the kernel to be in the middle of the array (at floor(N/2) + 1), so your t, for the purposes of building the kernel, should be t-floor(N/2) (this puts the 0 at the right location). Also, the kernel should be normalized to avoid changing the signal strength with the convolution. Just divide h by its sum.

2. But now, H = fft(h) will be wrong, because fft (and the DFT) expects the origin to be on the first element of the array. Use ifftshift to circularly shift the kernel array to put its origin on the first element: H = fft(ifftshift(h)).

3. For the manual DFT, you use the same matrix for the forward and the inverse transforms. You need to conjugate transpose the matrix to compute the inverse transform: y3 = (x * W .* H) * W'/N.

This is my corrected code (also some changes in the plotting to make it compatible with MATLAB):

N = 512;                  % number of points
t = 0:N-1;                % [0,1,2,...,N-1]
h = exp(-abs(t-floor(N/2))); % filter impulse response (this definition is symmetric around the origin, for a causal filter set the left half of h to zero)
h = h / sum(h);           % normalize
H = fft(ifftshift(h));    % filter frequency response
x = (1+t) .* sin(sqrt(1+t)); % (input signal of our choice)

y1 = conv(x,h,"same");    % Direct convolution
y2 = ifft(fft(x) .* H);   % FFT convolution

T  = transpose(t) * t;
W  = exp(2j*pi / N * T);  % DFT matrix
y3 = (x * W .* H) * W'/N; % "Manual" convolution

lw = 2;

clf
hold on
plot(t,x,        "displayname", "orig",                      "linewidth", lw+1)
plot(t,y1,       "displayname", "conv",   "linestyle", "--", "linewidth", lw)
plot(t,real(y2), "displayname", "fft",    "linestyle", ":",  "linewidth", lw)
plot(t,real(y3), "displayname", "manual", "linestyle", "-.", "linewidth", lw)
legend