3
$\begingroup$

If $x[n]$ is input signal, $y[m]$ is output signal, $h[m]$ is impulse response of the channel, then the following relation holds

$\lvert Y(\omega)\rvert = \lvert H(\omega)\rvert^2 . \lvert X(\omega)\rvert$

where $F(\omega)$ is the DFT of $f[n]$.

To verify this, I wrote following MATLAB code

x = randn(1,1000);
h = [1 2 3 4 5];
y = conv(x,h);
plot((abs(fft(h,1024)).^2).*(abs(fft(x,1024))));
hold on
plot(abs(fft(y,1024)),'r')

Result shows that second plot is at a much lower level than first plot. When I multiply second plot by $1024$ i.e. NFFT then both the plots overlap.

What is wrong with the original equation? I remember learning it that way in class.

enter image description here

$\endgroup$
  • $\begingroup$ Your claimed result $\lvert Y(\omega)\rvert = \lvert H(\omega)\rvert^2 . \lvert X(\omega)\rvert$ is incorrect and it is not to be wondered at that MATLAB does not confirm it. $\endgroup$ – Dilip Sarwate Jul 12 '14 at 11:12
4
$\begingroup$

Just change the code to the following:

x = randn(1,1000);
h = [1 2 3 4 5];
y = conv(x,h);
plot((abs(fft(h,1024))).*(abs(fft(x,1024)))); % It's |H(w)||X(w)|
hold on
plot(abs(fft(y,1024)),'--r')

By mistake you raised the DFT of an impulse response to the second power. You could see that magnitudes are bit off, but peaks and valleys are more-less at the same point.

enter image description here

Additionally I've noticed in my MATLAB 2014a, that you get an error after running your original code. That's because vector x was a row vector, and y a column one.


This part is to answer Drazick, claiming that this approach will produce some errors (already present in my solution).

In @jojek code there is no equality since he neglected the phase by applying the abs operator on both signals before the multiplication (As opposed to the Convolution Theorem).

If the filter had a more harsh phase his result would have much bigger error.

First of all please read my MATLAB code, maybe you will notice --r switch in second plot function. I intentionally used a dashed line, otherwise we would have just red plot and no blue at all.

But let's cut to the chase. Code is now slightly modified, as random signal is being filtered by long filter with random coefficients. This is what we get for comparison of $|H(\omega)X(\omega)|$ and $|H(\omega)||X(\omega)|$ with convolution in time domain:

enter image description here

Let's zoom in:

enter image description here

So all three are equal. Also script provides maximum difference between two methods to be around $10^{-13}$. Below script to reproduce the results:

clc, close all
x = randn(1,1000);
h = randn(1,200); % some crazy filter
y = conv(x,h);
N = length(y);

Y =  abs(fft(h,N)).*abs(fft(x,N)); % |H(w)||X(w)|
YY = abs(fft(h,N).*fft(x,N));      % |H(w)X(w)|

plot(YY,'go'); 
hold on
plot(Y,'b')    
plot(abs(fft(y,N)),'--r')
grid on
legend({'|H(w)X(w)|','|H(w)||X(w)|','conv'})

display(sprintf('Max error: %e', max(abs(YY-Y))))
$\endgroup$
  • $\begingroup$ But the Convolution theorem doesn't state the formula you implemented. You can't abs each of them by itself as you vanish the phase change due to the filtering. $\endgroup$ – Royi Jul 12 '14 at 7:13
  • $\begingroup$ see my updated code below. $\endgroup$ – Royi Jul 12 '14 at 7:34
  • 2
    $\begingroup$ @Drazick isnt abs(AB) = abs(A).abs(B).. see en.wikipedia.org/wiki/Absolute_value#equation_4 $\endgroup$ – user13107 Jul 12 '14 at 7:39
  • $\begingroup$ I mixed different case. Since in this case this is point wise multiplication it will work. Sorry. $\endgroup$ – Royi Jul 12 '14 at 10:01
4
$\begingroup$

I think you made a mistake between the power spectrum and the Fourier Transform. This is the right form of the Convolution Theory (Nothing is squared):

$$ \lvert Y(\omega)\rvert = \lvert H(\omega) X(\omega)\rvert $$

Try this and it will work for you on MATLAB.

MATLAB Code

% Calculating the Magnitude of the Filtered Signal

vFilterCoeff = [1; 2; 3; 4; 5];
vInputSignal = randn(1000, 1);

vOutputSginal = conv(vInputSignal, vFilterCoeff);

fftNumPoints = 1024;

vOutputSignalDftMagnitude = abs(fft(vOutputSginal, fftNumPoints));
vFilteredSignalDftMagnitude = abs(fft(vFilterCoeff, fftNumPoints) .* fft(vInputSignal, fftNumPoints));

hFigure = figure();
hAxes   = axes();
plot([vOutputSignalDftMagnitude, vFilteredSignalDftMagnitude]);

The result is:

enter image description here

As can be seen there's a full equality.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.