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

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

%input data, ]) {+ K  T7 `$ e& F9 O$ g
fa = 8000; %%signal frequency; @5 n6 z6 N0 g+ m) [: n
fs = 44100; %%44.1kHz sampling frequency
. Q' |3 J. d+ F* o, X, V3 ]n = 1:64;( {# y+ L$ ~4 C# l# N+ ~" \& T. d6 ^
x = sin(2*pi*fa*n/fs);1 P1 w. n% L8 S  u
lengthx = length(x);4 P% X* t. U$ c9 ]& U: G
t = n*1000/128/fs;- X# ?" D6 \3 L5 y
%plot(t,x);xlabel('time in ms');

% [y fs nbits] = wavread('acoustic.wav');2 j# P. ~3 _0 M' g) T* c/ [( R
% x = y';
" ^1 l( _% {& g6 Y% lx = length(x);
: Q4 P8 g9 h, k%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
) ^7 \7 _6 ~, e9 g%Resample parameters
: R8 m& n3 E4 E% VL1 = 8;
1 u% K0 Y: W+ G  M: C, S6 u0 SM1 = 7;4 E1 H5 U  Y. z4 e+ L' B9 x3 g
L2 = 4;
, A% A) ~  l, y& v: L- G/ u. E, {M2 = 3;! T, V& {+ j9 T. _3 v
L3 = 10;
" j4 m5 \) g- R7 eM3 = 7;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
6 R  J% @8 L8 _/ }/ Z%%filters
9 I% {2 _9 {* n5 Gload filter.mat4 @+ @' }: ?  {6 |! w, K: }
filter1 = b1;
" u9 {% b/ H" Y! d- @# k9 X, u4 K6 s$ sfilter2 = b2;8 x2 f' o; _9 M# z9 y& x
filter3 = b3;2 H6 A3 v" S1 U" W4 j" ]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%: i5 Q8 j) S. W( [* H; G" }$ ^
%Initialize the subfilters

%First Filter
' b" \  H1 S+ E& t' [upFilter1 = decompose(filter1,L1);5 H) I4 q) K/ N9 Q& \+ `" K
downFilter1 = decompose(filter1,M1);

%Second Filter
, J+ D. e- z( I3 @! z9 f0 RupFilter2 = decompose(filter2,L2);8 j7 U2 b/ h' Y9 w0 ^9 b
downFilter2 = decompose(filter2,M2);

%Second Filter
0 L- u- G( o( C1 M2 O4 P4 ^  QupFilter3 = decompose(filter3,L3);
" K& \( d6 f8 A" e4 D: T3 BdownFilter3 = decompose(filter3,M3);

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
' A1 Q+ J! Z' d%Polyphase Filtering

%First Filter
1 d( b. X. l4 Z7 |3 _) ^%
. N6 @, f, q! A1 D; X%%upsampling/ C6 {- R9 q. m) F& S0 L+ G
intermediatex1 = zeros(L1,length(x));
: [4 `) Y8 R' S5 Qx11 = zeros(L1,length(x)*L1);; U, u+ \5 X5 t2 Y
for i = 11
1 V/ A  \$ E: f7 z- ]8 Mintermediatex1(i, = filter(upFilter1(i,,1,x); %%do the filtering before upsampling: T1 V) z* h4 p. r
x11(i, = upsample(intermediatex1(i,:),L1); %%upsample by L, inserting L-1 zeros.5 b7 T7 L0 A2 U4 W! M
end

X1 = zeros(1,length(x)*L1);5 g' ^; u; t) k2 W" u
for i = 1:length(x)
% K9 n/ K& l6 S( v+ W9 a& dX1(L1*(i-1)+11*(i-1)+L1) = x11(:,(i-1)*L1+1)';2 ~7 v5 v' J& t) w2 H; z0 J( R# Z
end9 m% f8 g" \; X, c( X0 C
%%downsampling' Z/ ?- S0 h0 ]$ |2 }$ Y% H% ?7 c
X1_down = [X1 zeros(1,M1-mod(length(X1),M1))];
7 _( w' N6 J' \) N8 y6 I1 T& youtput1 = zeros(1,length(X1_down)/M1);% b7 e$ X# J. u
x1down = zeros(M1,length(X1_down)/M1);
9 [) o0 q( z/ zx11down = x1down;
4 A" q  t1 e( O, J0 T8 _: [7 ]for i = 1:M1
' V: Z% H# m6 k, [3 x$ N! d$ T7 vfor k = 1:length(X1_down)/M12 f: l( s* V0 ~3 b/ ^
x1down(i,k) = X1_down(M1*(k-1)+i);6 [: C3 ?" y$ R& E8 J4 ?) O
end
. E) i3 J" ?, _% sx11down(i,:) = filter(downFilter1(1,:),1,x1down(i,:));
; @: f2 W1 ?3 [' V; n- [output1 = output1 + x11down(i,:);" S2 z" ?  _8 f6 _4 d
end

%Second Filter
4 i+ _( c& s3 s%, `9 Y, j, F% Z! s
%%upsampling
3 l! t( r- q) g; Mintermediatex2 = zeros(L2,length(output1));( \4 B8 w& o2 M- p8 }2 D/ g% c
x22 = zeros(L2,length(output1)*L2);
! w& U8 n  }- e1 Xfor i = 12
. `. f9 _8 r$ D# p, p: J8 i! g8 fintermediatex2(i,:) = filter(upFilter2(i,:),1,output1); %%do the filtering before upsampling& o' L2 }; U& R% @  G
x22(i,:) = upsample(intermediatex2(i,:),L2); %%upsample by L, inserting L-1 zeros.' ^) `7 ]" \. d
end

X2 = zeros(1,length(output1)*L2);
, R9 n  d. J* v4 f: n. pfor i = 1:length(output1)
0 N# }: d  c+ |) F4 OX2(L2*(i-1)+1:L2*(i-1)+L2) = x22(:,(i-1)*L2+1)';
; D' @+ }( p$ g% Bend4 k4 y( z, j" K6 F1 z* f% O5 D5 W7 U
%%downsampling
/ G& F5 R, k9 R' MX2_down = [X2 zeros(1,M2-mod(length(X2),M2))];6 F3 u$ R$ Y% z  P
output2 = zeros(1,length(X2_down)/M2);
0 m6 U3 h$ Z' d  ux2down = zeros(M2,length(X2_down)/M2);3 e$ p+ F3 U. K3 \
x22down = x2down;
% N( J, P" ]- V, o3 J, G# xfor i = 1:M23 S* G6 e- C$ [' M* v9 k
for k = 1:length(X2_down)/M2
0 Y+ Q6 g- H# o8 j3 n- ^, tx2down(i,k) = X2_down(M2*(k-1)+i);
2 u( E& V& b9 K$ nend. m' F- o1 G4 m3 S1 w
x22down(i,:) = filter(downFilter2(1,:),1,x2down(i,:));
2 {& c/ I! x7 @$ @output2 = output2 + x22down(i,:);$ ?! R5 _+ v( r. z8 L9 y8 y% J, Z
end

%Third Filter
8 e* F4 u- }9 B5 J2 m& g+ M%: G6 n+ n& o3 Z  R4 g' j+ q  d
%%upsampling
- {+ n; |* f0 q1 D, F( `* }" M7 w0 gintermediatex3 = zeros(L3,length(output2));
, F/ \; H2 T2 |* w% ex33 = zeros(L3,length(output2)*L3);
; T* C- S2 D9 G% I% kfor i = 1:L3
- T' {2 t2 t+ w: @7 v9 m; V$ [! l: eintermediatex3(i,:) = filter(upFilter3(i,:),1,output2); %%do the filtering before upsampling8 f) x$ D1 R% V
x33(i,:) = upsample(intermediatex3(i,:),L3); %%upsample by L, inserting L-1 zeros.% W2 G0 h- J! s$ O
end

X3 = zeros(1,length(output2)*L3);
8 z1 R% [8 V1 a! f* _( jfor i = 1:length(output2)
1 j+ z( G. x# X" CX3(L3*(i-1)+1:L3*(i-1)+L3) = x33(:,(i-1)*L3+1)';
. X# M& f& D- {, d, M' z) Cend
9 O$ r* A7 @' m5 p2 ~0 [%%downsampling4 o1 S! J! v8 t4 y1 V% X
X3_down = [X3 zeros(1,M3-mod(length(X3),M3))];( H) v5 K$ ~  W- k  I6 K7 R
output3 = zeros(1,length(X3_down)/M3);
1 C$ n" y+ |- e- F; g: d. i% w, ]x3down = zeros(M3,length(X3_down)/M3);
/ l, O6 c+ @( i4 ]; e; Cx33down = x3down;* N2 L7 ?1 y& L* E: ^5 A. @+ _: o
for i = 1:M3
' A& W7 D9 k* a! D' S: Y9 }for k = 1:length(X3_down)/M36 t8 m# _* @4 ^2 y+ Z& Q
x3down(i,k) = X3_down(M3*(k-1)+i);
) t4 p) y0 m1 W6 oend
4 R7 J  c7 ?  e2 v1 @7 G5 W6 qx33down(i,:) = filter(downFilter3(1,:),1,x3down(i,:));4 R2 T4 q' R1 L- I9 o/ x
output3 = output3 + x33down(i,:);% R$ O" s0 o# {/ e3 w0 c0 O! Q- M5 w
end

output3 = output3/max(output3);* |$ w3 e+ G8 A3 g+ d, f
lout = length(output3);

nx = 1:length(output3);- J$ K: l% t8 X5 [1 d* u
xorigin = sin(2*pi*fa*nx/96000);
7 _. R3 X& f* G- z/ E# k' e" Uerror1 = xorigin - output3;7 S& V/ Z! W+ X
yin = [0 0 resample(x,320,147)];* f! O, M- x7 H8 L
error2 = xorigin - yin;

figure(1)
5 h5 q# D9 w4 b5 ?; y( Isubplot(3,1,1);plot(xorigin);title('True Signal');
. X( q9 K% p; F7 O$ X) hsubplot(3,1,2);plot(output3);title('Resampled Signal');
1 |) C. l. K8 `4 g( W. f) Psubplot(3,1,3);plot(error1);title('Error');

figure(2)- C' V0 C6 }' U/ {' l
subplot(3,1,1);plot(xorigin);title('True Signal');/ Y" G$ d3 b. P9 V8 A& S8 d2 r
subplot(3,1,2);plot(yin);title('Built-in Resampled Signal');
. `2 i1 s( M2 n/ Fsubplot(3,1,3);plot(error2);title('Error');

function x_decompose = decompose(x,factor)' t9 \0 Y- s* }  C: d" q: M8 A
lx = length(x);

if mod(lx,factor) ~= 0;
0 S3 [* Q+ d3 a' L" s( E, ?/ i5 sX = [x zeros(1,factor-mod(lx,factor))];5 I% j/ {5 \$ z' U* L) h
else# _3 K: F! T2 P  a& u. C0 _$ B1 v
X = x;
4 M6 T. G4 e0 U4 W1 i* Jend

x_decompose = zeros(factor,length(X)/factor);$ [& J3 q. r' M2 b, p9 H+ a# n
for i = 1:factor
2 l$ X- o0 d5 efor k = 1:length(X)/factor
% @$ K- |% W6 Z" _; vx_decompose(i,k) = X(factor*(k-1)+i);
$ A. r+ Q4 A1 ^# yend2 L* A. i  L3 E+ Y
end

end

%Filter Generation. ]& I7 P6 D2 e- n8 e& N/ R6 }
fs = 44100;
1 u+ J% c3 A9 o1 O& n. v. H. oL1 = 8;7 m) `% L  R; `9 U! y
M1 = 7;
( }2 g7 R3 a, y& k. oL2 = 4;& [0 ~; i( q. v% O' a+ T$ n  j" I5 F
M2 = 3;' u  M4 h( Q4 b! ]+ k
L3 = 10;
8 _& _( q1 c7 lM3 = 7;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%6 A9 y+ {# @" Q- M* l
%4 o0 E. a: C- M2 C5 F" h% h4 {5 w/ ]5 V
%%First Filter5 l# V9 h  Q; P
f1 = [20000 24000];
& i- f8 h* V* ~- x3 Ea1 = [1 0];: r7 j0 r' f+ S- @1 ?: g
dev1 = [0.01 10^-8];
3 G" @" g# q+ g" J. j- X! T/ zfs1 = fs*L1;  s  [, b7 h6 n; h8 \6 u0 f- c
[n1 fo1 ao1 w1] = firpmord(f1,a1,dev1,fs1);
* u% w# n2 X/ Q: K# M  Z% Rb1 = firpm(n1,fo1,ao1,w1);

%%Second Filter6 X( ?7 p/ l0 m
f2 = [20160 30240];
* f& |  R6 R2 H8 \* Ha2 = [1 0];% K5 ^. d4 N" B6 P& D
dev2 = [0.01 10^-8];# Y- C+ @2 E  T8 l  s, V8 g: N
fs2 = fs1*L2/M1;; G( h3 |: ?( X* _  ^- F: \6 a
[n2 fo2 ao2 w2] = firpmord(f2,a2,dev2,fs2);4 q$ K( {  T) Q0 M5 o' c
b2 = firpm(n2,fo2,ao2,w2);

%%Third Filter1 C& A8 n) `% c) T
f3 = [16800 50400];2 o& d1 W4 ]1 n9 I" S6 T' n2 q
a3 = [1 0];
) N, X6 J' k2 y- L' hdev3 = [0.01 10^-8];; v9 G. Q% N2 R7 s6 H! ~
fs3 = fs2*L3/M2;
7 Y1 z) Z8 `2 F' k[n3 fo3 ao3 w3] = firpmord(f3,a3,dev3,fs3);
; k, u5 [6 ?3 @; S+ f* t4 ob3 = firpm(n3,fo3,ao3,w3);