polyphase.m

信号多相滤波matlab 实现32 信道滤波

close all
clear all;
clc;

N=512;
M=32;
K=16;
n=1:N;

f1=600;
f2=5;
fs=2000;
t=-1:1/fs:1;
t1=0.2;
t2=0.5;
x1=sin(2*pi*f1*t);
% x2=square(2*pi*f2*t,20);
% x2=sawtooth(2*pi*f2*t,0.5);
x2=sin(2*pi*f2*t);
xx=x1+x2;
hh=firpm(N-1,[0 96 160 192*32]/(192*32),[1 1 0 0],[1 13]);
hh2=reshape(hh,M,K);

xk1=fft(x1);
xk2=fft(x2);
xxk=fft(xx);

figure(1)
subplot(2,3,1);plot(x1/max(abs(x1)));axis([0 2*fs -1 1]);xlabel('正弦波 x1');
subplot(2,3,2);plot(x2/max(abs(x2)));axis([0 2*fs -1 1]);xlabel('方波 x2');title('输入信号');
subplot(2,3,3);plot(xx/max(abs(xx)));axis([0 2*fs -1 1]);xlabel('xn=x1+x2');
subplot(2,3,4);plot(abs(xk1)/max(abs(xk1)));axis([0 2*fs 0 1]);
subplot(2,3,5);plot(abs(xk2)/max(abs(xk2)));axis([0 2*fs 0 1]);
subplot(2,3,6);plot(abs(xxk)/max(abs(xxk)));axis([0 2*fs 0 1]);
pause

figure(2)
subplot(3,1,1);
plot(xx/max(abs(xx)));hold on
plot(hh/max(abs(hh)),'r');
subplot(3,1,2);
ww=kaiser(1024,8)';
ww=ww/sum(ww);
plot((-0.5:1/1024:.5-1/1024)*32,fftshift(20*log10(abs(fft(hh,1024)))),'r');axis([-16 16 -100 0]);hold on
plot((-0.5:1/1024:.5-1/1024)*32,fftshift(20*log10(abs(fft(xx(1:1024).*ww,1024)))));  
subplot(3,1,3);
plot((-0.5:1/1024:.5-1/1024)*32,fftshift(20*log10(abs(fft(xx(1:1024),1024)))));axis([-16 16 0 50]);
pause

Q=N*4
gg=reshape(hh,M,K); % 抽取
reg=zeros(M,K);
n2=1;
for nn=1:M:Q
reg(:,2:K)=reg(:,1:K-1); %右移
reg(:,1)=flipud(xx(nn:nn+M-1)');  %排序
for mm=1:M
vv(mm)=reg(mm,:)*gg(mm,:)';
end
yy(:,n2)=fft(abs(vv))';
n2=n2+1;
end

p=0;
figure(3)
subplot(5,2,1)
plot(abs(real(yy(p+1,:)))/max(abs(real(yy(p+1,:)))));title('Time Series, Channels 1-5');
subplot(5,2,3)
plot(abs(real(yy(p+2,:)))/max(abs(real(yy(p+2,:)))));%grid;axis([0 32 -1 1]);
subplot(5,2,5)
plot(abs(real(yy(p+3,:)))/max(abs(real(yy(p+3,:)))));%grid;axis([0 32 -1 1]);
subplot(5,2,7)
plot(abs(real(yy(p+4,:)))/max(abs(real(yy(p+4,:)))));%grid;axis([0 32 -1 1]);
subplot(5,2,9)
plot(abs(real(yy(p+5,:)))/max(abs(real(yy(p+5,:)))));%grid;axis([0 32 -1 1]);

% ww=kaiser(200,8)';
% ww=ww/sum(ww);
subplot(5,2,2)
fft(yy(p+1,:));
plot((-0.5:1/512:.5-1/512)*32,(20*log10(abs(fft(yy(p+1,:),512)))));
subplot(5,2,4)
plot((-0.5:1/512:.5-1/512)*32,(20*log10(abs(fft(yy(p+2,:),512)))));
subplot(5,2,6)
plot((-0.5:1/512:.5-1/512)*32,(20*log10(abs(fft(yy(p+3,:),512)))));
subplot(5,2,8)
plot((-0.5:1/512:.5-1/512)*32,(20*log10(abs(fft(yy(p+4,:),512)))));
subplot(5,2,10)
plot((-0.5:1/512:.5-1/512)*32,(20*log10(abs(fft(yy(p+5,:),512)))));