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

loveeatmore

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

%input data
fa = 8000; %%signal frequency
0 X; ?5 ^1 X/ Y  c/ {+ q* H% ?fs = 44100; %%44.1kHz sampling frequency
5 r$ k  Q3 ?& f/ O' L6 J) B8 Jn = 1:64;
/ |6 b& |2 l# y" Q" @x = sin(2*pi*fa*n/fs);
lengthx = length(x);
+ T7 q, z0 R  ^8 ]+ q% c( bt = n*1000/128/fs;
, G* A1 M- \) |2 o%plot(t,x);xlabel('time in ms');

% [y fs nbits] = wavread('acoustic.wav');
+ D( C2 Q, D: B/ z& D3 ~( z3 n! l% x = y';
+ u7 V: X: _9 \( F. I4 W9 \# o% lx = length(x);
0 S" r9 d: ~; m8 Q%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Resample parameters
L1 = 8;
7 @+ l7 ^' ~: U$ eM1 = 7;
  S; j/ V: @& S% G: r8 P5 kL2 = 4;
M2 = 3;
L3 = 10;
- k- ^% k$ p' XM3 = 7;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
/ A: a6 F1 i( p3 n# h%%filters
load filter.mat
filter1 = b1;
filter2 = b2;
+ X/ a, G- |; F6 zfilter3 = b3;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Initialize the subfilters

%First Filter
upFilter1 = decompose(filter1,L1);
downFilter1 = decompose(filter1,M1);

%Second Filter
4 a5 `7 r. @8 \4 ^upFilter2 = decompose(filter2,L2);
downFilter2 = decompose(filter2,M2);

%Second Filter
$ g' G* @+ h0 |4 tupFilter3 = decompose(filter3,L3);
downFilter3 = decompose(filter3,M3);

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Polyphase Filtering

%First Filter
; Z5 Q2 B0 Y5 T- @%
/ H4 Z- u( e/ {' _  Z%%upsampling
intermediatex1 = zeros(L1,length(x));
x11 = zeros(L1,length(x)*L1);
7 Q7 t: w) q: B2 V/ Pfor i = 11
intermediatex1(i, = filter(upFilter1(i,,1,x); %%do the filtering before upsampling
x11(i, = upsample(intermediatex1(i,:),L1); %%upsample by L, inserting L-1 zeros.
3 B) ~2 L+ E: uend

X1 = zeros(1,length(x)*L1);
for i = 1:length(x)
X1(L1*(i-1)+11*(i-1)+L1) = x11(:,(i-1)*L1+1)';
/ h& m6 Z, D% y; H  V& ^; J& Cend
5 s1 s, b0 [: t  \5 o/ i3 y2 Y- y- M%%downsampling
X1_down = [X1 zeros(1,M1-mod(length(X1),M1))];
output1 = zeros(1,length(X1_down)/M1);
. H! _% [) O$ {: e7 H2 O  s/ t7 Rx1down = zeros(M1,length(X1_down)/M1);
- a/ m$ D, P  e3 \: l  Lx11down = x1down;
" R/ b0 v! U  o! Z; Dfor i = 1:M1
+ Y1 l7 I, Q" T- s  gfor k = 1:length(X1_down)/M1
6 r3 m, h2 M$ }6 a- [. px1down(i,k) = X1_down(M1*(k-1)+i);
& p! _- \% s6 s9 H) s+ ^3 `; H5 Rend
3 X6 N' o2 ?: w2 h* B8 Wx11down(i,:) = filter(downFilter1(1,:),1,x1down(i,:));
. A( L- P1 x! w8 o# b; ooutput1 = output1 + x11down(i,:);
end

%Second Filter
%
%%upsampling
intermediatex2 = zeros(L2,length(output1));
x22 = zeros(L2,length(output1)*L2);
" b# F5 K+ u- M7 }1 `# w/ Tfor i = 12
4 N, C) _( e! j2 ^/ D4 ~) Tintermediatex2(i,:) = filter(upFilter2(i,:),1,output1); %%do the filtering before upsampling
x22(i,:) = upsample(intermediatex2(i,:),L2); %%upsample by L, inserting L-1 zeros.
end

X2 = zeros(1,length(output1)*L2);
+ _0 x3 v. a# N% tfor i = 1:length(output1)
' k$ P7 i3 X0 F$ GX2(L2*(i-1)+1:L2*(i-1)+L2) = x22(:,(i-1)*L2+1)';
end
%%downsampling
( u7 `) u& a+ |0 ?) ~1 fX2_down = [X2 zeros(1,M2-mod(length(X2),M2))];
0 |! j. k8 i9 E3 ]output2 = zeros(1,length(X2_down)/M2);
x2down = zeros(M2,length(X2_down)/M2);
/ j4 y. L9 _& ax22down = x2down;
for i = 1:M2
$ z' K' G  m) A/ u; N' a/ ufor k = 1:length(X2_down)/M2
& K# u! a+ W- _' K" `8 w4 K4 Ux2down(i,k) = X2_down(M2*(k-1)+i);
2 T/ ~2 M' @- I+ ]2 F' ~end
' R) S8 }6 I3 U. {5 _+ \x22down(i,:) = filter(downFilter2(1,:),1,x2down(i,:));
0 K% ?* S6 U* Zoutput2 = output2 + x22down(i,:);
end

%Third Filter
%
* ]4 `% X  t+ [8 [9 O%%upsampling
2 N. n: ^1 D" ]9 K) N. fintermediatex3 = zeros(L3,length(output2));
' Q) H6 c/ I: p& s  Yx33 = zeros(L3,length(output2)*L3);
" O9 j3 d2 y4 c9 gfor i = 1:L3
intermediatex3(i,:) = filter(upFilter3(i,:),1,output2); %%do the filtering before upsampling
) [% z) n/ N: P! d# {x33(i,:) = upsample(intermediatex3(i,:),L3); %%upsample by L, inserting L-1 zeros.
end

X3 = zeros(1,length(output2)*L3);
for i = 1:length(output2)
, r  W/ y) E, D1 P/ yX3(L3*(i-1)+1:L3*(i-1)+L3) = x33(:,(i-1)*L3+1)';
# D  V/ S# r" I  Q; [1 kend
5 k' X2 h3 S" U$ p8 Z%%downsampling
  b7 N6 b! j9 T0 MX3_down = [X3 zeros(1,M3-mod(length(X3),M3))];
output3 = zeros(1,length(X3_down)/M3);
" `$ s; }6 }+ p! Y7 T% h. E9 ?x3down = zeros(M3,length(X3_down)/M3);
6 _3 T2 i0 ^: j- E: }x33down = x3down;
5 U! I2 w9 w- L. tfor i = 1:M3
, t4 S4 f: r' `9 Bfor k = 1:length(X3_down)/M3
x3down(i,k) = X3_down(M3*(k-1)+i);
end
& P6 v3 G; w3 `( t  E/ w1 j$ c) yx33down(i,:) = filter(downFilter3(1,:),1,x3down(i,:));
, r- z6 ?  m0 s1 _/ Poutput3 = output3 + x33down(i,:);
end

output3 = output3/max(output3);
6 B) @1 T) w# e5 Z' Tlout = length(output3);

nx = 1:length(output3);
xorigin = sin(2*pi*fa*nx/96000);
error1 = xorigin - output3;
yin = [0 0 resample(x,320,147)];
error2 = xorigin - yin;

figure(1)
subplot(3,1,1);plot(xorigin);title('True Signal');
subplot(3,1,2);plot(output3);title('Resampled Signal');
7 Q  o# h4 V- J" y* a3 msubplot(3,1,3);plot(error1);title('Error');

figure(2)
subplot(3,1,1);plot(xorigin);title('True Signal');
* X1 {3 q# d5 c! }5 vsubplot(3,1,2);plot(yin);title('Built-in Resampled Signal');
  h' |* `$ s! d1 ^! r- \subplot(3,1,3);plot(error2);title('Error');

---

function x_decompose = decompose(x,factor)
lx = length(x);

if mod(lx,factor) ~= 0;
X = [x zeros(1,factor-mod(lx,factor))];
* \) u: t# I# |, y0 ~  Q, C5 @else
) Q1 K( U* w) M3 n, e7 @X = x;
( ~  X2 U$ D+ l/ U: Hend

x_decompose = zeros(factor,length(X)/factor);
for i = 1:factor
$ ~& g1 U, T( }; i' ofor k = 1:length(X)/factor
3 }  l& V3 I+ S3 ex_decompose(i,k) = X(factor*(k-1)+i);
end
end

end

%Filter Generation
0 z1 j0 f+ A; ~0 G  Yfs = 44100;
L1 = 8;
M1 = 7;
L2 = 4;
+ Y3 ^1 W& e  E# CM2 = 3;
L3 = 10;
5 g5 \1 }" K. @! g: z" ^! \* p/ F& RM3 = 7;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
- a* U5 F6 h9 d%%First Filter
; z( k, p- G9 T$ M7 ~% O& ]6 I/ a% Of1 = [20000 24000];
a1 = [1 0];
dev1 = [0.01 10^-8];
fs1 = fs*L1;
[n1 fo1 ao1 w1] = firpmord(f1,a1,dev1,fs1);
b1 = firpm(n1,fo1,ao1,w1);

%%Second Filter
f2 = [20160 30240];
# w) W! f7 c' _  Oa2 = [1 0];
dev2 = [0.01 10^-8];
) l) {" W. ~8 ^fs2 = fs1*L2/M1;
- j4 p7 T5 L" u. s' B  V$ _6 `/ k[n2 fo2 ao2 w2] = firpmord(f2,a2,dev2,fs2);
b2 = firpm(n2,fo2,ao2,w2);

%%Third Filter
( x: A& b5 \, e9 h7 @6 Cf3 = [16800 50400];
a3 = [1 0];
dev3 = [0.01 10^-8];
7 V' d! V4 C2 d" X1 ffs3 = fs2*L3/M2;
! ?  I4 \. A5 }# ]8 j  p[n3 fo3 ao3 w3] = firpmord(f3,a3,dev3,fs3);
b3 = firpm(n3,fo3,ao3,w3);

木棉花_MM

这个……