2015-08-18 69 views
0

我正在嘗試在Matlab中建立一個帶有不同頻段的濾波器。過濾後,低頻被推遲,因爲你可以在下面的圖片看到:濾波後發生低頻移位

過濾前: Before filtering

過濾後: After filtering

你能幫助我理解爲什麼會這樣?

我使用的代碼如下:

[x,Fs] = audioread('drum-loop.wav'); 

% define constants - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

ntaps = 255;                % number of filter-taps 

% sampling frequencies for each branch 
fs_A = 48000;     % 
fs_B = 24000;     % 
fs_C = 6000;     % 
fs_D = 1500;     % 

% decimation for each branch 
dec_A = 1;      % 
dec_B = 2;      % 
dec_C = 4;      % 
dec_D = 4;      % 

% decimation filters - - - - - - - - - - - - - - - - - - - - - - - - - - - 
% ref: https://es.mathworks.com/help/dsp/ref/mfilt.html 
decfilter_B = mfilt.firdecim(dec_B);          % decimate signal for branch B 
x_B = filter(decfilter_B, x);            % apply decimation filter 

decfilter_C = mfilt.firdecim(dec_C);          % decimate signal for branch C 
x_C = filter(decfilter_C, x_B);           % apply decimation filter 

decfilter_D = mfilt.firdecim(dec_D);          % decimate signal for branch D 
x_D = filter(decfilter_D, x_C);           % apply decimation filter 

% interpolation filters - - - - - - - - - - - - - - - - - - - - - - - - - 
% ref: https://es.mathworks.com/help/dsp/ref/mfilt.html 
intfilter_B = mfilt.firinterp(dec_B);          % decimate signal for branch B 
intfilter_C = mfilt.firinterp(dec_C);          % decimate signal for branch C 
intfilter_D = mfilt.firinterp(dec_D);          % decimate signal for branch D 

% define filters - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
% refs: 
% https://es.mathworks.com/help/signal/ref/fir1.html 
% http://es.mathworks.com/matlabcentral/answers/17260-fir1-basic-question 
% http://es.mathworks.com/help/pdf_doc/signal/signal_tb.pdf 

ovlp1 = 1000; 
ovlp2 = 500; 
ovlp3 = 250; 

% branch A: band pass filter. 
coeff_A = [fs_B/2-ovlp1 (fs_A/2-1)]/(fs_A/2); 
bp_A = fir1(ntaps,coeff_A); 

% branch B: band pass filter. 
coeff_B = [fs_C/2-ovlp2 (fs_B/2-1)]/(fs_B/2); 
bp_B = fir1(ntaps,coeff_B); 

% branch C: band pass filter. 
coeff_C = [fs_D/2-ovlp3 (fs_C/2-1)]/(fs_C/2); 
bp_C = fir1(ntaps,coeff_C); 

% branch D: band pass filter. 
coeff_D = (fs_D/2-1)/(fs_D/2); 
bp_D = fir1(ntaps,coeff_D,'low'); 

% apply filters 
x_A_f = filter(bp_A,1,x); 
x_B_f = filter(bp_B,1,x_B); 
x_C_f = filter(bp_C,1,x_C); 
x_D_f = filter(bp_D,1,x_D); 

% summation filter 
y = x_A_f + filter(intfilter_B, x_B_f) + ... 
    filter(intfilter_B, filter(intfilter_C, x_C_f)) + ... 
    filter(intfilter_B, filter(intfilter_C, filter(intfilter_D, x_D_f))); 
y = y/max(y);  % normalize 

% % compute fft 
figure; 
subplot(3,1,1); 
plot(x); 
NFFT = 2^nextpow2(length(x));           % Next power of 2 from length of y (in samples) 
Y = fft(x,NFFT)/Fs;  
f = Fs/2*linspace(0,1,NFFT/2+1); 
subplot(3,1,2);spectrogram(x,blackman(128),60,128,1e3,'yaxis') 
subplot(3,1,3);plot(f,2*angle(Y(1:NFFT/2+1))) 

% % compute fft 
figure; 
subplot(3,1,1); 
plot(y); 
NFFT = 2^nextpow2(length(y));           % Next power of 2 from length of y (in samples) 
Y = fft(y,NFFT)/Fs;  
f = Fs/2*linspace(0,1,NFFT/2+1); 
subplot(3,1,2);spectrogram(y,blackman(128),60,128,1e3,'yaxis') 
subplot(3,1,3);plot(f,2*angle(Y(1:NFFT/2+1))) 

audiowrite('output.wav',y,48000); 

回答

1

您的過濾器引入的相移。使用filtfilt而不是filter