MATLAB: Calculate the frequency range and plot the unwrapped phase spectrum

frequency rangeparseval's theoremunwrapped phase spectrum

I am given a DTFT, how can I calculate the frequency range where 65% of the total energy is contained? I did unwrap the phase, but I don't know if it's right(because my unwrapped plot looks the same as original). Can anybody help? I will appreciate it. Here it is my script:
% Read in the numerator and denominator coefficients                                      
num = [0.1323 0.1323*0.1444 -0.1323*0.4519 0.1323*0.1444 0.1323];
den = [1 0.1386 0.8258 0.1393 0.4153];
% Compute the frequency response
w = 0:0.01:2*pi;
h = freqz(num, den, w);
% Plot the frequency response
subplot(2,2,1)
plot(w/pi,real(h));grid
title('Real part')
xlabel('\omega/\pi'); ylabel('Amplitude')
subplot(2,2,2)
plot(w/pi,imag(h));grid
title('Imaginary part')
xlabel('\omega/\pi'); ylabel('Amplitude')
subplot(2,2,3)
plot(w/pi,abs(h));grid
title('Magnitude Spectrum')
xlabel('\omega/\pi'); ylabel('Magnitude')
subplot(2,2,4)
plot(w/pi,angle(h));grid
title('Phase Spectrum')
xlabel('\omega/\pi'); ylabel('Phase, radians')
%unwrapping the phase
unwrap=unwrap(angle(h),0.8);
figure, subplot(211)
plot(w/pi,angle(h));grid
title('Phase spectrum');
xlabel('\omega/\pi'); ylabel('Phase, radians');
subplot(212);
plot(w/pi,unwrap);grid
title('Unwrapped phase spectrum');            
xlabel('\omega/\pi'); ylabel('Phase, radians');

Best Answer

Honglei makes two good points in his comments.
I'll assume you are talking about the distribution of the squared-magnitudes of the frequency response. If you obtain the frequency response of your filter, you can normalize that frequency response (it's just a vector) so that the sum of the magnitude-squared frequency responses equal 1.
   [H W]=freqz(num,den);
   H = H./norm(H);
   pwr = abs(H).^2;
   % sum(pwr) is 1
   sum(pwr(W>=0 & W<=pi))
Note in the last line I have used logical indexing on the frequency vector. If you just plot the magnitude-squared frequency response
 plot(W,pwr)
You clearly see that essentially all the power is between approx. 1 and 2.5 radians. You can confirm that with
 sum(pwr(f>=1 & f<=2.5))
The problem with the way that you have phrased your question is that you can choose many such intervals that have "65% of the total energy". For example:
 sum(pwr(W>=0 & W<=1.72))
will give you something close to 0.65, but so will:
sum(pwr(W>=1.4 & W<=pi))
So which interval do you choose?
Related Question