Why does freqz function make a phase difference?

I was dealing with my DSP project that this issue suddenly crossed my way.

I was using two different method to calculated the Fourier transform of a vector.

Then figured out that there is a "phase difference" between the results, gained by "Freqz" and "equation".

equation: ; n varies from zer to infinity.

It seems like that freqz has a delay that makes the results different.

But when you insert a 0 as the first element of the vector; the result would be the same as the one calculated by "equation" without changing the vector.

Can anyone tell me that what is the reason behind this?

clc;clear;close all;
%calcuating fourier transform of x by freqz function.
A=1;
x=[1 2 3 4];
[H,W1]=freqz(A,x,1000);
subplot(211);
plot(W1,(angle(H)),'b-');
hold on;
%calcuating fourier transform of x by the equation.
W=0:0.001:1*pi;
N=length(x);
xw=zeros(1,length(W));
    for n=1:N
        xw=xw+x(n).*exp(1i.*n.*W);
    end
hold on;
plot(W,(angle(xw)),'r-'); grid;
xlabel('W'); ylabel('phase');
legend('angle calculated by freqz','angle calculated by the equation');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%calcuating fourier transform of x by freqz function.
x=[1 2 3 4];  x2=[0 1 2 3 4];
[H,W2]=freqz(A,x2,1000);
subplot(212);
plot(W2,(angle(H)),'b-');
hold on;
%calcuating fourier transform of x by the equation by x NOT x2.
W=0:0.001:1*pi;
N=length(x);
xw=zeros(1,length(W));
    for n=1:N
        xw=xw+x(n).*exp(1i.*n.*W);
    end
plot(W,(angle(xw)),'r-'); grid;
xlabel('W'); ylabel('phase');
legend('angle calculated by freqz','angle calculated by the equation');

Answer 1

[Paul]

Hi Mohammad,

See code below for corrections

%calcuating fourier transform of x by freqz function.
A=1;
x=[1 2 3 4];

The arguments to freqz need to be reversed to compute the DTFT of x as in the equation in the question.

%[H,W1]=freqz(A,x,1000);
[H,W1]=freqz(x,A,1000);

subplot(211);
plot(W1,(angle(H)),'b-');
hold on;
%calcuating fourier transform of x by the equation.
W=0:0.001:1*pi;
N=length(x);
xw=zeros(1,length(W));
    for n=1:N

As in the equation in the question, the argument to exp() needs a negative sign. And as stated in the question, the n in the sum starts at 0, so we have to subtract 1 from the n used in the for loop that has to start from 1 because Matlab uses 1-based indexing. Or you could make the for loop go from 0 to N-1 and then use x(n+1).

%        xw=xw+x(n).*exp(1i.*n.*W);
        xw=xw+x(n).*exp(-1i.*(n-1).*W);
    end
hold on;
plot(W,(angle(xw)),'r-'); grid;
xlabel('W'); ylabel('phase');
legend('angle calculated by freqz','angle calculated by the equation');

[Mohammad]

Hi Paul

first of all thanks for your comment.

I corrected the power of exp and added a negative sign.

And for [H,W]=freqz(B,A,N) I put the numerator vector B= x and denominator equal to 1.

But as in your opinion you mentioned that the delay is due to Matlab 1-base indexing; if l use a 0-base indexing programming language like Python (by assuming that it has such a freqz function); l wouldn't have such a problem? And this is all due to Matlab structure?

[Paul]

With the corrected code, I don't know what you are refererring to as "the delay" and there is no problem as far as I can tell.

I was assuming that your finite duration signal was mathematically defined as:

x[n] = n+1 for 0 <= n < = 3

x[n] = 0, otherwise.

i.e., x[0] = 1, x[1] = 2, x[2] = 3, x[3] = 4

The Discrete Time Fourier Transform of x[n] is, by your equation,

x[0]*exp(-1j*0*w) + x[1]*exp(-1j*1*w) + x[2]*exp(-1j*2*w) + x[3]*exp(-1j*3*w) (1)

Now we define the Matlab variable X as (I've capitalized it to make it easier to distinguish)

X = [1 2 3 4]

The for loop summation written in Matlab code computes (1), but has to account for the fact that Matlab uses 1-based indexing into Matlab vectors and so the relationship between the mathemtical object x[n] and the Matlab vector X is

X(n) = x[n-1]

If you want to see that explicitly, we can "unroll" four iterations of the loop to compute the DTFT

xw = xw + X(1).*exp(-1i.*(0).*W); % n = 1

xw = xw + X(2).*exp(-1i.*(1).*W); % n = 2

xw = xw + X(3).*exp(-1i.*(2).*W); % n = 3

xw = xw + X(4).*exp(-1i.*(3).*W); % n = 4

If you're using a language that allows zero-based indexing, then the loop could be written as:

for n=0:N-1 % assuming N is set to 4

xw = xw+X(n).*exp(-1i.*n.*W);

end

and you should get the same result.

This call

freqz(X,1,1000)

computes (1) at 1000 angular frequencies around the upper half of the unit circle as can be seen here, where we have to keep in mind that the Matlab programming variables X and b use one-based indexing, hence,

X(i) = b(i) = x[i-1]

So there is no delay (I'm still not sure what you meant by that); there is only the need to translate mathematical equations into the Matlab programming language under the constraints of the language and using the tools that the language offers. Which would be true of any language.