关于重采样的问题，麻烦大神帮忙看看

自己在matlab里面编的一个重采样程序，但是效果很差，而且输出前一段数据有错，麻烦路过的大神帮我看看？

%input data( D( L" g- C6 g9 s6 p
fa = 8000; %%signal frequency
& d/ c! R, W* pfs = 44100; %%44.1kHz sampling frequency. z) C7 {2 [6 S
n = 1:64;& E- g0 H$ V$ Z0 F9 A5 n) R
x = sin(2*pi*fa*n/fs);) q/ [% `+ u; `9 a1 `
lengthx = length(x);
, D+ ], ^, w# ]9 Y" `0 xt = n*1000/128/fs;$ E/ ^  r9 s  F/ e' T- p# [
%plot(t,x);xlabel('time in ms');

% [y fs nbits] = wavread('acoustic.wav');
" {& A1 J; V2 s# @+ f% x = y';
1 _. G) N+ F% Q* v$ G7 z% lx = length(x);) n8 K$ N0 L; d0 ]0 F8 |. G' k
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
; q  n4 |/ J! k' y6 F2 D%Resample parameters0 _( N5 ~9 }% T: F) ?! x0 a
L1 = 8;( o6 N: N7 i4 J2 s/ E0 b: K
M1 = 7;9 p8 @% H$ }+ h
L2 = 4;/ Q% p5 N3 [5 H6 n7 E6 t
M2 = 3;1 v+ Y) z8 p8 L8 X; O
L3 = 10;
0 \- _- v8 ?$ v9 h' @M3 = 7;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
7 j0 y5 }* ]' _/ c# }5 p%%filters4 m8 I+ A7 F$ g% e2 ^
load filter.mat1 Z( I9 @5 q: [( W4 i4 _
filter1 = b1;
) Z. ?* o- x: U. D0 ofilter2 = b2;
# P# B( V: T, E0 K3 S: Rfilter3 = b3;! M2 ]5 z# ~) L+ P9 _# R8 W; z+ z! X
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%& l: ~( ]! @/ F8 f' o
%Initialize the subfilters

%First Filter# r, m: x, C% S$ N* Z# t5 b
upFilter1 = decompose(filter1,L1);4 e6 ]8 |) J' b" w0 x( M' V" R7 ^9 y2 p
downFilter1 = decompose(filter1,M1);

%Second Filter
7 J; q3 K. m, H' UupFilter2 = decompose(filter2,L2);
9 ]' m7 p, n. C& ?5 L0 ?' VdownFilter2 = decompose(filter2,M2);

%Second Filter
/ {% `$ p8 s) ~- SupFilter3 = decompose(filter3,L3);  v8 s. C% N* y
downFilter3 = decompose(filter3,M3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 t/ K5 X: t( Z5 p0 H: u%Polyphase Filtering

%First Filter
6 @8 i- y! b8 N%
+ l4 Q; g" e1 l  d0 m( Y* x%%upsampling; \; b8 y0 p* n: o: @. A" O* V% ]
intermediatex1 = zeros(L1,length(x));1 g+ ?. N; u$ Q+ ]
x11 = zeros(L1,length(x)*L1);, b7 P( Z1 j( F7 ]8 c* g1 d
for i = 11
6 f, t3 K/ m0 |! Y# {intermediatex1(i, = filter(upFilter1(i,,1,x); %%do the filtering before upsampling, _1 V5 `/ M' M
x11(i, = upsample(intermediatex1(i,:),L1); %%upsample by L, inserting L-1 zeros.1 n5 |% u* d* A* X5 Z; a
end

X1 = zeros(1,length(x)*L1);
3 D6 e/ y* S  T6 r; T0 `for i = 1:length(x)- \  {8 ]& ?! G9 r' r5 W
X1(L1*(i-1)+11*(i-1)+L1) = x11(:,(i-1)*L1+1)';
6 L+ u! D+ Q8 `- r. [6 x* N" Y/ Uend
7 t3 v( Y) Q5 u! B%%downsampling: P; @- R% s8 v+ ]; \: G! Z7 {
X1_down = [X1 zeros(1,M1-mod(length(X1),M1))];
) _8 \" _3 W0 Y. foutput1 = zeros(1,length(X1_down)/M1);
! N% a) ~; j" W" Kx1down = zeros(M1,length(X1_down)/M1);
. i, a+ D- K# G1 j0 L6 yx11down = x1down;, n7 ]7 X( E! V* a
for i = 1:M1
$ P$ ~- w! X* \4 l3 O* Ffor k = 1:length(X1_down)/M17 c: q' D5 X+ G% N. \& t' T: o
x1down(i,k) = X1_down(M1*(k-1)+i);
6 v: |. F' E  G2 ~% }' l$ _$ Dend2 [/ u  l- l$ R( z
x11down(i,:) = filter(downFilter1(1,:),1,x1down(i,:));0 {! k2 i9 J9 Y- a5 o) K! }9 v3 G2 ?
output1 = output1 + x11down(i,:);' S: i1 B; R( O5 \
end

%Second Filter$ a1 {, Q, Q0 ]5 w
%) K# P3 s7 W! K# N$ t
%%upsampling5 s7 D) s. g0 k1 ]
intermediatex2 = zeros(L2,length(output1));
+ n/ |  B& W, m' _- kx22 = zeros(L2,length(output1)*L2);
, ~9 p# i* Q, w+ h  k! E3 K/ ]for i = 123 A+ Y5 X, x; c8 [) n; a6 {
intermediatex2(i,:) = filter(upFilter2(i,:),1,output1); %%do the filtering before upsampling
' v# A# D$ ], F7 _x22(i,:) = upsample(intermediatex2(i,:),L2); %%upsample by L, inserting L-1 zeros.- e) s) Z/ k+ I: l. _2 `
end

X2 = zeros(1,length(output1)*L2);% m- X8 R5 O# }2 M$ t
for i = 1:length(output1)
8 y* q8 l' g; k6 N( pX2(L2*(i-1)+1:L2*(i-1)+L2) = x22(:,(i-1)*L2+1)';& K2 S+ k# K: w7 X# M3 G& o
end* l. d7 ?( u- S
%%downsampling; R. P6 c8 P; t; o* O
X2_down = [X2 zeros(1,M2-mod(length(X2),M2))];
& W$ i% M5 \" I( h+ y$ uoutput2 = zeros(1,length(X2_down)/M2);
2 }2 l7 J2 z4 O1 s, s, y+ Ux2down = zeros(M2,length(X2_down)/M2);$ W8 k. T' |6 {. D
x22down = x2down;
& k: W) d" f2 Nfor i = 1:M2# U3 ^. Y. ~7 M( A- u
for k = 1:length(X2_down)/M2
' ^2 h/ d9 S! Ix2down(i,k) = X2_down(M2*(k-1)+i);  n; u1 X# S0 ~6 S0 F/ k$ ?
end8 ?9 _9 ]: W  D; N. v+ ^9 ]0 a
x22down(i,:) = filter(downFilter2(1,:),1,x2down(i,:));4 W. N; G" L6 n3 F0 J" C
output2 = output2 + x22down(i,:);- [! X6 u! h# ^7 W& `
end

%Third Filter
! @$ i$ L$ d# `- S0 \# m+ d+ }%
8 S" [" ]3 y$ J5 ]" {3 \- Q+ h; U%%upsampling
5 O* f2 j& ~  Yintermediatex3 = zeros(L3,length(output2));. j% ~) C- Q( B9 i4 ~: A4 n- I! d
x33 = zeros(L3,length(output2)*L3);
+ [1 V9 ?1 o6 O. `$ Y/ y$ i! O/ |for i = 1:L3# O' I/ R. d4 H" [. C9 s
intermediatex3(i,:) = filter(upFilter3(i,:),1,output2); %%do the filtering before upsampling
1 h4 l6 l. ?3 R9 Hx33(i,:) = upsample(intermediatex3(i,:),L3); %%upsample by L, inserting L-1 zeros.& y$ ^; p! W4 G' B
end

X3 = zeros(1,length(output2)*L3);
6 s  E& e+ @5 Z% Nfor i = 1:length(output2)6 B9 W/ u$ Q6 j+ |& i6 e0 b( |  |
X3(L3*(i-1)+1:L3*(i-1)+L3) = x33(:,(i-1)*L3+1)';" ^, g  G# S& G# m& q% S
end
( B3 F) x# J- B) E% w%%downsampling- K8 _6 q; w/ k
X3_down = [X3 zeros(1,M3-mod(length(X3),M3))];
0 h2 W& @% a, `& I& G4 u# B, Soutput3 = zeros(1,length(X3_down)/M3);
# N2 ]. }, z/ @+ N/ Q7 l: y& @  [x3down = zeros(M3,length(X3_down)/M3);
" E, i, l3 J2 p/ C" ~x33down = x3down;
3 k8 C) A8 V. c3 w7 b& Yfor i = 1:M3" x4 u( n( K- f5 X* ^  M) g$ N, L( Y
for k = 1:length(X3_down)/M3
, I+ N1 c5 l9 C* I' }( yx3down(i,k) = X3_down(M3*(k-1)+i);1 Y2 ~; e! l" v& k+ ^' d6 \% l$ m% P
end' t& U+ V) m/ R/ a7 @
x33down(i,:) = filter(downFilter3(1,:),1,x3down(i,:));7 L( S2 F3 S( Z1 p% b/ N6 I
output3 = output3 + x33down(i,:);# G5 `/ b% p5 T7 T8 }
end

output3 = output3/max(output3);
/ n1 U& }# e# j  H6 }( G" |lout = length(output3);

nx = 1:length(output3);
  d, i" a0 v" i! c- zxorigin = sin(2*pi*fa*nx/96000);
7 \  _' m9 S0 u0 f- s$ c# kerror1 = xorigin - output3;
2 U8 C$ ]8 R. y, d& P% e& _yin = [0 0 resample(x,320,147)];4 N( F8 a; P" C! V, m& x
error2 = xorigin - yin;

figure(1)5 {3 j0 i& i  P5 o
subplot(3,1,1);plot(xorigin);title('True Signal');5 O! h4 E$ E; ?* O0 n
subplot(3,1,2);plot(output3);title('Resampled Signal');
( l7 e6 M( l' F& Gsubplot(3,1,3);plot(error1);title('Error');

figure(2)# ?8 c* Z/ |) _" G, n- r
subplot(3,1,1);plot(xorigin);title('True Signal');
6 n9 e* f( ]# M/ h+ ksubplot(3,1,2);plot(yin);title('Built-in Resampled Signal');
! ]& F; s! z$ w; U7 ^9 D0 h; P+ }" Esubplot(3,1,3);plot(error2);title('Error');

function x_decompose = decompose(x,factor)0 P% z3 Y. R$ m- ^8 C9 I* C4 _
lx = length(x);

if mod(lx,factor) ~= 0;
! j9 m/ l* u5 R! V4 hX = [x zeros(1,factor-mod(lx,factor))];1 }4 k$ x9 ]- f2 X& L
else
$ q  h* C$ c2 `X = x;
: _; y$ E- g; U" F: v* W; c: fend

x_decompose = zeros(factor,length(X)/factor);. @; n) J5 u' j+ u! o4 Q
for i = 1:factor6 ?* o5 Z) H2 d& A: e( {
for k = 1:length(X)/factor
: U' n; K9 u: k$ X  Px_decompose(i,k) = X(factor*(k-1)+i);  z% `1 q6 [/ n- V! Q9 c" F. W
end
1 _+ Y' x* \( N- e! d( }end

end

%Filter Generation
" ~' I$ @- l: `5 g( \, I) V+ Afs = 44100;8 |9 Z7 s" b- r* l$ p
L1 = 8;
( c( H; J* A1 a" f' xM1 = 7;
: H3 a, n$ O  n1 YL2 = 4;$ [- G% @$ g2 S# F
M2 = 3;0 e+ K/ J" O: U5 E1 _- A
L3 = 10;' q7 K' L3 f2 Z& b9 Q+ J
M3 = 7;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%: {" A5 I7 E. M( u/ j3 Q
%. e8 q& N# l; ~3 u3 Q
%%First Filter: R- y4 t; K) m5 V' @
f1 = [20000 24000];+ ?) Q) K# }8 C0 O, ^5 G' `
a1 = [1 0];" F! U( w2 i2 s* m4 ?  ^, p
dev1 = [0.01 10^-8];
( o; g! G1 Y5 ^7 q$ Mfs1 = fs*L1;9 Y5 M2 [8 T1 H$ }2 f/ O
[n1 fo1 ao1 w1] = firpmord(f1,a1,dev1,fs1);+ g9 |1 L7 `0 S1 I* {1 A
b1 = firpm(n1,fo1,ao1,w1);

%%Second Filter
7 C. p" n2 L5 j* [1 G/ tf2 = [20160 30240];
3 e: D' d4 L  ta2 = [1 0];9 V# Z9 s. X7 W) [0 o4 M$ J
dev2 = [0.01 10^-8];
. g7 P& z+ M/ j. n, E# Xfs2 = fs1*L2/M1;' I( B9 V4 B1 Q! }
[n2 fo2 ao2 w2] = firpmord(f2,a2,dev2,fs2);$ M& A: |, Y4 V! c  q, \
b2 = firpm(n2,fo2,ao2,w2);

%%Third Filter1 c* I" G  u! W* A$ i+ P' h' D
f3 = [16800 50400];9 M8 E4 h$ i: U. ]+ |
a3 = [1 0];& j2 S  j( f6 J% ?+ z
dev3 = [0.01 10^-8];
% T8 L- H& [- ofs3 = fs2*L3/M2;' f: J: @# |6 E. F' c  F. C
[n3 fo3 ao3 w3] = firpmord(f3,a3,dev3,fs3);
" M9 H, u  y; w. g, Jb3 = firpm(n3,fo3,ao3,w3);