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

loveeatmore

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

%input data
fa = 8000; %%signal frequency
) J0 S3 q) a0 Zfs = 44100; %%44.1kHz sampling frequency
  g  J8 D& Z+ U; j3 [9 b# wn = 1:64;
x = sin(2*pi*fa*n/fs);
+ @3 ^  s5 n8 @lengthx = length(x);
t = n*1000/128/fs;
9 Z4 P+ T- {4 y- I+ W: ?1 ?, f( J%plot(t,x);xlabel('time in ms');

% [y fs nbits] = wavread('acoustic.wav');
# \& [+ C: N6 J. ]2 O  z% x = y';
' N% F, y% ?* E( E* z8 Y/ o$ Q% lx = length(x);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
9 i& P$ @2 Y* @1 S3 S* Z4 z( F( ]. x%Resample parameters
L1 = 8;
M1 = 7;
: h7 t1 K6 n+ `1 U* P& }, y; QL2 = 4;
M2 = 3;
L3 = 10;
7 f/ k8 c& A7 r' ?M3 = 7;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%filters
load filter.mat
& i; q/ K7 X/ w* p2 b3 Ffilter1 = b1;
: J$ S, n5 v4 u0 Mfilter2 = b2;
1 i% }$ x$ o; g; yfilter3 = b3;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
7 s% ?$ ~9 o' B9 B7 i%Initialize the subfilters

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

%Second Filter
5 V' s6 g, X/ u5 mupFilter2 = decompose(filter2,L2);
downFilter2 = decompose(filter2,M2);

%Second Filter
6 d) l$ N0 T8 }& n; hupFilter3 = decompose(filter3,L3);
downFilter3 = decompose(filter3,M3);

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- }4 |8 g( Y5 Y%Polyphase Filtering

%First Filter
! C% {- b" q5 e%
7 p" H2 G7 `4 k0 @& u%%upsampling
intermediatex1 = zeros(L1,length(x));
x11 = zeros(L1,length(x)*L1);
for i = 11
4 Q, f) h# E2 q" T' Hintermediatex1(i, = filter(upFilter1(i,,1,x); %%do the filtering before upsampling
) x' u, ]+ I# ^5 F+ d) M* bx11(i, = upsample(intermediatex1(i,:),L1); %%upsample by L, inserting L-1 zeros.
end

X1 = zeros(1,length(x)*L1);
for i = 1:length(x)
, T2 T6 O' g' A/ t% w" @X1(L1*(i-1)+11*(i-1)+L1) = x11(:,(i-1)*L1+1)';
- A; ?2 @0 l: s5 f( C8 P( fend
+ V, e1 i3 q" |( v  i* K* b9 ]%%downsampling
X1_down = [X1 zeros(1,M1-mod(length(X1),M1))];
output1 = zeros(1,length(X1_down)/M1);
6 c7 Q" r( Q/ L- q- Cx1down = zeros(M1,length(X1_down)/M1);
x11down = x1down;
for i = 1:M1
: C" T' h& R% Yfor k = 1:length(X1_down)/M1
x1down(i,k) = X1_down(M1*(k-1)+i);
6 Y; P* Y5 y8 D" v: _/ f; M9 pend
x11down(i,:) = filter(downFilter1(1,:),1,x1down(i,:));
output1 = output1 + x11down(i,:);
0 a' N- L3 i5 `9 P3 E" s" b7 L" fend

%Second Filter
%
%%upsampling
intermediatex2 = zeros(L2,length(output1));
x22 = zeros(L2,length(output1)*L2);
for i = 12
intermediatex2(i,:) = filter(upFilter2(i,:),1,output1); %%do the filtering before upsampling
x22(i,:) = upsample(intermediatex2(i,:),L2); %%upsample by L, inserting L-1 zeros.
: F% S& j( j7 B3 eend

X2 = zeros(1,length(output1)*L2);
: k: y8 v: y9 W& o9 B  Lfor i = 1:length(output1)
1 B6 z7 {0 H- T0 c- H& jX2(L2*(i-1)+1:L2*(i-1)+L2) = x22(:,(i-1)*L2+1)';
5 ?* R+ E) [$ K3 oend
%%downsampling
X2_down = [X2 zeros(1,M2-mod(length(X2),M2))];
4 M- X3 L( t' X* _' t' ?% Noutput2 = zeros(1,length(X2_down)/M2);
( A4 \4 x% C3 s9 a4 D- G9 P3 i' h$ jx2down = zeros(M2,length(X2_down)/M2);
x22down = x2down;
for i = 1:M2
$ z1 P. Z1 z8 X6 jfor k = 1:length(X2_down)/M2
x2down(i,k) = X2_down(M2*(k-1)+i);
end
x22down(i,:) = filter(downFilter2(1,:),1,x2down(i,:));
# q7 n- P/ |* loutput2 = output2 + x22down(i,:);
4 d4 j0 n1 k; j0 F* t* jend

%Third Filter
7 Y$ R+ z; ?4 Q2 E%
1 f" M: e* q9 z7 h! h( m: p3 v%%upsampling
+ H# Y& \8 r# I: J  Rintermediatex3 = zeros(L3,length(output2));
x33 = zeros(L3,length(output2)*L3);
for i = 1:L3
intermediatex3(i,:) = filter(upFilter3(i,:),1,output2); %%do the filtering before upsampling
# s1 d& A* V6 q# w4 e# s% b, Sx33(i,:) = upsample(intermediatex3(i,:),L3); %%upsample by L, inserting L-1 zeros.
( P, F! h9 h/ X, ~9 S# `$ Lend

X3 = zeros(1,length(output2)*L3);
for i = 1:length(output2)
5 K( x, N1 n% m3 x0 a1 p8 ^X3(L3*(i-1)+1:L3*(i-1)+L3) = x33(:,(i-1)*L3+1)';
: J, }9 O9 N' C9 d. x) M3 iend
%%downsampling
3 k7 {5 _; J7 P6 f! PX3_down = [X3 zeros(1,M3-mod(length(X3),M3))];
& ^: Z+ g$ f: M, L) j* koutput3 = zeros(1,length(X3_down)/M3);
5 ~, E: a8 {% N- Rx3down = zeros(M3,length(X3_down)/M3);
4 [) y2 v% a1 Q8 w" A/ M7 gx33down = x3down;
for i = 1:M3
for k = 1:length(X3_down)/M3
x3down(i,k) = X3_down(M3*(k-1)+i);
end
x33down(i,:) = filter(downFilter3(1,:),1,x3down(i,:));
: }4 z) K. h. a6 y4 M' goutput3 = output3 + x33down(i,:);
end

output3 = output3/max(output3);
lout = length(output3);

nx = 1:length(output3);
6 y% d5 [) P9 |- b' dxorigin = sin(2*pi*fa*nx/96000);
. f6 f. G- ^, W0 uerror1 = xorigin - output3;
* H( B: _. f5 T& Iyin = [0 0 resample(x,320,147)];
% D7 m. F, y+ P8 T' o1 Derror2 = xorigin - yin;

figure(1)
7 q7 u( W" Z5 Y8 Xsubplot(3,1,1);plot(xorigin);title('True Signal');
subplot(3,1,2);plot(output3);title('Resampled Signal');
subplot(3,1,3);plot(error1);title('Error');

figure(2)
7 H; B9 G; E1 _: b7 usubplot(3,1,1);plot(xorigin);title('True Signal');
) @' U9 V0 Z7 `. H; Wsubplot(3,1,2);plot(yin);title('Built-in Resampled Signal');
- M, w, C8 s  c5 y$ rsubplot(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))];
else
0 l+ c$ i* D( n, S. I( B6 K  D& {& G5 oX = x;
9 q9 ?  P( v5 Qend

x_decompose = zeros(factor,length(X)/factor);
( J, H) d3 Q- A" s# W7 ?! cfor i = 1:factor
for k = 1:length(X)/factor
' p5 g8 u" J2 K( F6 o0 s* dx_decompose(i,k) = X(factor*(k-1)+i);
end
) I- k1 [# V' j2 B9 L' M( Gend

end

%Filter Generation
: I$ X1 [0 N0 k6 m0 G: I1 O  K9 ]fs = 44100;
# |4 M' |8 N  N, H; y" T" X( u1 tL1 = 8;
4 s1 i) z6 F1 w9 X) eM1 = 7;
4 u, ?6 o0 Z- K8 W" sL2 = 4;
M2 = 3;
L3 = 10;
M3 = 7;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4 w4 r, }5 q! _1 j4 b! y' ?$ I$ L0 V4 \%
4 Z1 W( K2 f! m/ q%%First Filter
f1 = [20000 24000];
: h  j  ^' V7 E( `( m  ~5 a( s1 za1 = [1 0];
dev1 = [0.01 10^-8];
fs1 = fs*L1;
[n1 fo1 ao1 w1] = firpmord(f1,a1,dev1,fs1);
" |- Y9 c. C6 {% G6 T- z( r" w' ub1 = firpm(n1,fo1,ao1,w1);

%%Second Filter
f2 = [20160 30240];
a2 = [1 0];
9 b9 h. i- B$ t% udev2 = [0.01 10^-8];
fs2 = fs1*L2/M1;
[n2 fo2 ao2 w2] = firpmord(f2,a2,dev2,fs2);
b2 = firpm(n2,fo2,ao2,w2);

%%Third Filter
f3 = [16800 50400];
a3 = [1 0];
% b* a0 h2 i6 v" C$ E$ K, Pdev3 = [0.01 10^-8];
; U+ A; P; _1 @4 D; sfs3 = fs2*L3/M2;
[n3 fo3 ao3 w3] = firpmord(f3,a3,dev3,fs3);
b3 = firpm(n3,fo3,ao3,w3);

8 M7 c* K. X# z* Z8 u. G2 o
; P' R& @" ^9 J, Y

木棉花_MM

这个……