filterbank8.m

8 通道双正交余弦调制滤波器组实现

%************************************************************%
% 8 通道双正交余弦调制滤波器组实现
% 原型滤波器的实现 利用matlab的滤波器设计工具
% 滤波器的采样频率为 Fs:20000Hz 设计的低通原型的通带频率为625Hz 阻带衰减80db
% 滤波器的系数长度为512
clear; clc;
fid=fopen('fcoff.txt','r'); 
if(fid==-1)
    fprintf('cannot open the file!\n');
    return;
end
h0_str=fscanf(fid,'%s');
fclose(fid);
h0=str2num(h0_str);
clear h0_str;

% 8个分析滤波器和8个综合滤波器系数实现
Fs=10000;
M=8;                                 % 通道数
Delay = 2*2*M;                       % Delay = 2*s*M+d,其中s为正整数,0<= d <= 2*M-1
Theta = pi/4;
n = length(h0);
Hkn = zeros(M,n);
Gkn = zeros(M,n);
for k = 1:M
    for i = 1:n
        Hkn(k,i)=2*h0(i)*cos(((k-1)+0.5)*(i-1-Delay/2)*pi/M + (-1)^(k-1)*pi/4);
        Gkn(k,i)=2*h0(i)*cos(((k-1)+0.5)*(i-1-Delay/2)*pi/M - (-1)^(k-1)*pi/4);
    end
end

%************************************
% 时域系数观察
figure(1)
subplot(311)
stem(h0)
title('原型滤波器系数h0')
subplot(312)
stem(Hkn(1,:));
title('分析滤波器的系数H0')
subplot(313)
stem(Gkn(1,:));
title('综合滤波器的系数G0')
%***********************************
% 功率互补验证
sumh = zeros(1,512);
for k=1:M
    [H(k,:),w] = freqz(Hkn(k,:),1,512,Fs);
    h(k,:) = 20*log10(abs(H(k,:)+eps));
    sumh = abs(H(k,:)).^2 + sumh;
end
figure(2)
subplot(211)
for k=1:M
    plot(w/Fs,h(k,:));
    hold on;
end
hold off;
title('分析滤波器组的幅频响应')
sum=10*log10(sumh);
subplot(212);
plot(w/Fs,sum);grid
title('分析滤波器幅频响应和')
%**********************************************
% 测试信号生成
t = 0:1/Fs:0.3;
x=exp(j*2*pi*500*t)+0.1*randn(size(t))+exp(j*2*pi*800*t)+exp(j*2*pi*3000*t)+exp(j*2*pi*4000*t);
%x=exp(j*2*pi*500*t)+0.1*randn(size(t))+exp(j*2*pi*800*t)+exp(j*2*pi*300*t)+exp(j*2*pi*600*t);
%**********************************************
% 信号通过滤波器组
for k=1:M
    Hkdeci = mfilt.firdecim(M,Hkn(k,:));
    xk(k,:) = filter(Hkdeci,x);
    Xkinter =mfilt.firinterp(M,Gkn(k,:));
    xg(k,:) = filter(Xkinter,xk(k,:));
end
x_end = zeros(1,length(xg));                    % 恢复出来的信号
for k=1:M
    x_end = xg(k,:)+x_end;
end
figure(3)
subplot(211);
plot((1:1000),x(1:1000));
title('输入信号')
subplot(212);
plot((1:999),x_end(1:999))
title('输出信号')


figure(4)
subplot(211);
h=spectrum.periodogram;               
psd(h,x);                                        %输入信号频谱图
title('输入信号频谱')
subplot(212);
h=spectrum.periodogram;
psd(h,x_end);                                     %输出信号频谱图
title('输出信号频谱')