信息处理 - 通过线性相位滤波器测量群延迟 - 吾爱随笔录

我很难解决相位/群延迟问题，它与您开始采样的位置相对应，因此我建立了一个系统来测试其中一些问题。系统将 2 个正弦波通过具有线性相位的低通 FIR 滤波器，然后我从输出数据的某个任意点开始取一个 NFFT 长度 fft，然后测量每个正弦波的相位。信号流程如下：

[正弦波] => [FIR（长度 65）] => [FFT] => [相位测量]

我正在尝试使用可以从 fft 中的相位信息计算组延迟的事实：

但这就是我迷路的地方，因为您在输出端开始采样的位置似乎会影响相位测量。即你不能从任意点开始。但是在一个黑盒系统中，您无法确定从什么正确的输出样本开始，您如何才能准确地测量群延迟？下面是我的脚本。我的群延迟应该是 32，但是根据变量out_start_idx的群延迟变化很大。不应该是不管你从哪里开始对输出进行采样，两个音之间的相位变化应该仍然是恒定的吗？

close all
clear all

fs = 1e6;

%GROUP DELAY of filter (N-1)/2 = 32
%since h = length 65
h = remez(2^6, [0 fs*.05 fs*.15 fs/2]/(fs/2), [1 1 0 0]); 

%arbitrarily large fft
NFFT = 2^19; 

bin_space  = fs/NFFT;

num_tones = 2; 
space = 1000; % spacing in bins
for nn = 1:num_tones
 freqs(nn) = space*(nn)*bin_space;
end

%random starting phases between -pi/2 pi/2
phi_start = rand(1,num_tones)*pi-pi/2; 
t = 0:999999; 

sig = 0; 
for nn = 1:num_tones
 sig = sig + exp(1i*2*pi*freqs(nn)*t/fs + phi_start(nn));
end

out = filter(h,1, [sig zeros(1,length(h)+ 10)]); 

%{
figure
subplot(211)
plot((-.5:1/NFFT:.5-1/NFFT)*fs, 20*log10(fftshift(abs(fft(h,NFFT)))))
hold on
plot((-.5:1/NFFT:.5-1/NFFT)*fs, 20*log10(fftshift(abs(fft(sig,NFFT)))))
hold off
subplot(212)
plot((-.5:1/NFFT:.5-1/NFFT)*fs, 20*log10(fftshift(abs(fft(out,NFFT)))))
hold on
plot((-.5:1/NFFT:.5-1/NFFT)*fs, 20*log10(fftshift(abs(fft(h,NFFT)))))
hold off
%}
f = 0:fs/NFFT:fs-fs/NFFT;

% take NFFT amount of samples from arbitrary starting point 
% in output
out_start_idx = 1; 
data_test = out(out_start_idx:out_start_idx+NFFT-1); 

% get index of tones
idx_tones = (space:space:space*(num_tones)) + 1 

% take fft of data
fft_out = fft(data_test,NFFT);

% get complex output from each tone bin
tone_IQ = fft_out(idx_tones); 

% calculate phase of each tone
phi_tones = angle(tone_IQ)*180/pi;

% calculate group delay
phi_2 = diff(fliplr(phi_tones))
phi_1 = diff(fliplr(phi_start)*180/pi); 
tau = (-1/360) * (phi_2 - phi_1)./diff(freqs)
grp_delay = tau*fs

% maybe do some extra processing
phi_out = phi_tones; 

% linear fit calculated phases
p = polyfit(freqs,phi_out,1);
y = polyval(p, freqs); 

% calculate deviation from linear phase
deviation_from_linear_phase = phi_out -y; 

% plot deviation / etc
figure
plot(freqs, y)
hold on
plot(freqs, phi_out, 'r*')

figure
stem(deviation_from_linear_phase)