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

loveeatmore

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

%input data
4 u! M- a1 j) A7 kfa = 8000; %%signal frequency
fs = 44100; %%44.1kHz sampling frequency
n = 1:64;
1 {) d1 o3 u' Cx = sin(2*pi*fa*n/fs);
1 N  v! [! F/ l- @% G2 q( S% olengthx = length(x);
t = n*1000/128/fs;
9 i' I) X& I6 G& H* [- w%plot(t,x);xlabel('time in ms');

% [y fs nbits] = wavread('acoustic.wav');
4 }& ~' @" ^0 o6 u; `) N% x = y';
+ `/ E, a+ M( l5 \- g: l0 V3 I+ Q( D% lx = length(x);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
; N, X4 g& Z0 f: w7 b# J" \! _3 O: t( G%Resample parameters
- ]7 a6 u2 s2 N% |! ~L1 = 8;
' w7 J0 v: T6 U7 S, O% [M1 = 7;
! l( t: @  C) a& G& j" iL2 = 4;
M2 = 3;
& ^7 G- R2 i0 j- @L3 = 10;
0 e5 l; s# N2 G  uM3 = 7;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
* j' ^2 G  Q8 }9 z$ l3 l- _% P%%filters
# K6 s2 m1 e% r/ Z- M8 F/ Zload filter.mat
filter1 = b1;
. D1 |, [+ X: x% d- a+ Pfilter2 = b2;
filter3 = b3;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Initialize the subfilters

%First Filter
1 Y; u) {) ~5 q5 uupFilter1 = decompose(filter1,L1);
1 C' Q. V# n0 v  ^downFilter1 = decompose(filter1,M1);

%Second Filter
upFilter2 = decompose(filter2,L2);
/ M( j; `; i6 }. ]: `/ {, QdownFilter2 = decompose(filter2,M2);

%Second Filter
upFilter3 = decompose(filter3,L3);
downFilter3 = decompose(filter3,M3);

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
* D5 y$ b; n* P/ t! \  |%Polyphase Filtering

%First Filter
# t6 W6 h+ K% p2 h7 |% h9 \%
' p7 W* i0 j. n" h2 u$ e! ^3 d1 O%%upsampling
intermediatex1 = zeros(L1,length(x));
x11 = zeros(L1,length(x)*L1);
+ _1 G9 k4 m+ I' b% dfor i = 11
4 z! u8 S; N0 C* `: a6 ?! I) W4 H7 |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.
2 h2 |! m1 H3 t, r+ i7 xend

X1 = zeros(1,length(x)*L1);
for i = 1:length(x)
X1(L1*(i-1)+11*(i-1)+L1) = x11(:,(i-1)*L1+1)';
end
%%downsampling
4 [% \- ?9 l4 p1 \X1_down = [X1 zeros(1,M1-mod(length(X1),M1))];
" U9 x6 k3 B4 routput1 = zeros(1,length(X1_down)/M1);
. p! D: U0 ?1 ^+ E! f5 ~x1down = zeros(M1,length(X1_down)/M1);
x11down = x1down;
9 F0 {2 h5 m2 l* W( W* P# Tfor i = 1:M1
for k = 1:length(X1_down)/M1
' x! N# O. R3 p4 \9 kx1down(i,k) = X1_down(M1*(k-1)+i);
$ Z! i; o+ q3 K3 |% G8 rend
: g4 \, X; P7 G7 o$ e2 Yx11down(i,:) = filter(downFilter1(1,:),1,x1down(i,:));
output1 = output1 + x11down(i,:);
end

%Second Filter
5 C) L8 [, A. G6 b" y  g( x; Y7 }%
' w* G8 P! ]) `, ^5 N8 l1 A) e%%upsampling
intermediatex2 = zeros(L2,length(output1));
x22 = zeros(L2,length(output1)*L2);
, B& O( b: @( J+ ]2 ]( ?+ }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.
end

X2 = zeros(1,length(output1)*L2);
+ T0 [9 v6 }6 g% Q" x2 F9 _for i = 1:length(output1)
, Y% ?0 S& r6 `: \. nX2(L2*(i-1)+1:L2*(i-1)+L2) = x22(:,(i-1)*L2+1)';
) S- B9 {$ }, m7 ^8 b. t- Nend
%%downsampling
X2_down = [X2 zeros(1,M2-mod(length(X2),M2))];
  c. g0 v+ u& \output2 = zeros(1,length(X2_down)/M2);
x2down = zeros(M2,length(X2_down)/M2);
x22down = x2down;
7 _/ k7 g2 `) w' `7 Tfor i = 1:M2
for k = 1:length(X2_down)/M2
5 q% l, L0 J* y/ P" t& Ix2down(i,k) = X2_down(M2*(k-1)+i);
# I& A; I1 y( j1 ~0 send
$ D1 w; D! Z! L1 M! {/ x- zx22down(i,:) = filter(downFilter2(1,:),1,x2down(i,:));
4 [5 J& W% G, b1 B& routput2 = output2 + x22down(i,:);
end

%Third Filter
4 ?/ k7 `* p! u# a5 i7 C%
) _# b6 H  J, T& Z1 c%%upsampling
* G. O  y+ v; |& g6 L+ f' |9 f% l: Vintermediatex3 = 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
- R4 z7 ~  f1 g- v( P- f6 @# Zx33(i,:) = upsample(intermediatex3(i,:),L3); %%upsample by L, inserting L-1 zeros.
end

X3 = zeros(1,length(output2)*L3);
for i = 1:length(output2)
X3(L3*(i-1)+1:L3*(i-1)+L3) = x33(:,(i-1)*L3+1)';
" H5 f) t9 Q4 T" U7 B1 R+ d$ b' N7 hend
%%downsampling
X3_down = [X3 zeros(1,M3-mod(length(X3),M3))];
output3 = zeros(1,length(X3_down)/M3);
x3down = zeros(M3,length(X3_down)/M3);
x33down = x3down;
2 ^& h' G: z5 hfor i = 1:M3
0 [$ J; H" c8 r3 ^for k = 1:length(X3_down)/M3
x3down(i,k) = X3_down(M3*(k-1)+i);
end
+ W0 k' ]2 K- S, `$ ?) sx33down(i,:) = filter(downFilter3(1,:),1,x3down(i,:));
output3 = output3 + x33down(i,:);
% {0 t" _; p" y& q# K+ U6 h+ [end

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

nx = 1:length(output3);
xorigin = sin(2*pi*fa*nx/96000);
' V4 b5 d3 W, {" u+ X% P% z4 Ferror1 = xorigin - output3;
yin = [0 0 resample(x,320,147)];
error2 = xorigin - yin;

figure(1)
subplot(3,1,1);plot(xorigin);title('True Signal');
6 [# }. @& ]0 \0 r  ~- ysubplot(3,1,2);plot(output3);title('Resampled Signal');
subplot(3,1,3);plot(error1);title('Error');

figure(2)
subplot(3,1,1);plot(xorigin);title('True Signal');
subplot(3,1,2);plot(yin);title('Built-in Resampled Signal');
subplot(3,1,3);plot(error2);title('Error');

---

function x_decompose = decompose(x,factor)
' y9 P- U+ |; q; L- U, O1 ulx = length(x);

if mod(lx,factor) ~= 0;
X = [x zeros(1,factor-mod(lx,factor))];
; ]* u6 k/ R, Q6 Selse
X = x;
- r7 l2 ~  a* u" G/ h# j) bend

x_decompose = zeros(factor,length(X)/factor);
) ~# M+ m5 F" ?* G9 Afor i = 1:factor
+ l: J( h; A8 L/ O( Tfor k = 1:length(X)/factor
x_decompose(i,k) = X(factor*(k-1)+i);
7 ^, K- I, T0 nend
/ U$ o$ @: ?# n0 ]; bend

end

%Filter Generation
fs = 44100;
- g1 w" V9 x/ D4 wL1 = 8;
& j; u+ X9 G- O3 j- oM1 = 7;
L2 = 4;
M2 = 3;
# f8 m% G+ j: Q& ^L3 = 10;
M3 = 7;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
" @2 q0 ?1 U$ A+ ?  x%
%%First Filter
f1 = [20000 24000];
a1 = [1 0];
0 a- M/ R, Q  O4 e  f2 kdev1 = [0.01 10^-8];
fs1 = fs*L1;
) X7 Q  g0 R: C+ M5 X[n1 fo1 ao1 w1] = firpmord(f1,a1,dev1,fs1);
2 D, l/ f* ~" |. sb1 = firpm(n1,fo1,ao1,w1);

%%Second Filter
( E1 x0 v3 G( _+ ~" g$ g3 Qf2 = [20160 30240];
a2 = [1 0];
4 L7 p& b6 c9 F$ U: t' \dev2 = [0.01 10^-8];
7 i" s, b" H& q2 q) P4 Sfs2 = fs1*L2/M1;
* ^6 Q" B1 v/ Y; A) R[n2 fo2 ao2 w2] = firpmord(f2,a2,dev2,fs2);
b2 = firpm(n2,fo2,ao2,w2);

%%Third Filter
f3 = [16800 50400];
' {3 @3 Q) K9 Z2 b( s0 @a3 = [1 0];
dev3 = [0.01 10^-8];
fs3 = fs2*L3/M2;
+ z. t4 P1 m+ J; c[n3 fo3 ao3 w3] = firpmord(f3,a3,dev3,fs3);
  I; ]9 O, _; o4 F2 Ob3 = firpm(n3,fo3,ao3,w3);

  G2 ?/ c* `8 w. U

木棉花_MM

这个……