Manual Extraction of Filter Phase Response with MATLAB using sweep sine

Question

I'm currently at a loss for why my code seems to not be working. I'm doing a sine sweep on a filter to calculate its freq and phase response. The frequency response checks out, however the phase response doesn't look right when compared to MATLABs freqz function. Plus since I know the filter is a comb, I can tell its wrong as well.

Here's what I've got so far.

%   Filter Coefficients
b = [1];
a = [1, 0, 0, 0, 0, 0, 0, 0, 0, -0.96];

%  sweep is my sine wave that goes from 0->fs/2
filtered_sweep = filter(b, a, sweep);

% This section correctly Grabs and plots the Freq Response
fft_sweep = fft(filtered_sweep);
freq_response = abs(fft_sweep);
freq_response = 20*log10(freq_response/1);      % Convert to dB scale
freq_response = freq_response(1:end/2);         % Cut to pi fs/2    
figure;
plot(freq_response);

And Now here is where I try and get the Phase two different ways with the same wrong result. From what I understand to get the phase you just need the angle of the complex values. In MATLAB I thought angle() should work, but just to be safe I tried to manually compute the angle myself using atan with the imag and real components.

% Pre Allocate Vector
% Manual Angle Extraction without using Angle()
phase = zeros(1, length(fft_sweep));
for i=1:length(fft_sweep)
    phase(i) = atan2(imag(fft_sweep(i)), real(fft_sweep(i)));
end

% Using angle()
phase_response = angle(fft_sweep);
phase_response = phase_response(1:end/2);
figure;
plot(unwrap(phase_response));
grid on;

Both methods give the same incorrect phase response. Any help would be appreciated.

Answer 1

You missing one, very important step. In order to get the impulse response of your filter while using the sweep sine, you must get the output of your system $h[n]$ to excitation with sweep-sine signal $x[n]$:

$y[n]=h[n]\star x[n]$

In your case that's the filtering with filter function. Now in order to get the impulse response you must convolve your signal with the inverse filter $f[n]$. Generally that's the time-inverted sweep signal (although special kind of pre-conditioning, amplitude scaling, etc. might be done as well):

$h[n]=y[n]\star f[n]$

Now you have your impulse response that you can use for further analysis. Please take a look at this answer tackling very similar problem.

For more theory please see:

A. Farina - Advancements in Impulse Response Measurements by Sine Sweeps

A. Farina - Simultaneous Measurement of Impulse Response and Distortion with a Swept-Sine Technique