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

loveeatmore

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

%input data
fa = 8000; %%signal frequency
fs = 44100; %%44.1kHz sampling frequency
" u9 G2 ]9 d1 k/ Y! a, P8 ?n = 1:64;
6 y$ u, V( \( v* v  Z. L. Gx = sin(2*pi*fa*n/fs);
lengthx = length(x);
t = n*1000/128/fs;
%plot(t,x);xlabel('time in ms');

% [y fs nbits] = wavread('acoustic.wav');
, D) o( a- v4 o0 ~- ~9 e% x = y';
% lx = length(x);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! u/ Q  i. Y( o3 o: ]+ B%Resample parameters
L1 = 8;
M1 = 7;
% f- L3 j  {6 X; M9 T1 O' |L2 = 4;
9 x: ~8 W  b5 n: UM2 = 3;
L3 = 10;
7 C) R/ u$ o9 D" k. Q  xM3 = 7;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
" z6 N. T. ]- h/ a: L$ J, a%%filters
load filter.mat
- X" n7 t& X+ T" C( C+ Ofilter1 = b1;
$ D! ~' I. Z! A- ofilter2 = b2;
) s* e% y+ g( i( Tfilter3 = b3;
5 H" u: t9 T% g+ g. M( I) {: [$ q%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Initialize the subfilters

%First Filter
upFilter1 = decompose(filter1,L1);
( t8 u& v& c6 V+ ]downFilter1 = decompose(filter1,M1);

%Second Filter
upFilter2 = decompose(filter2,L2);
; U& c7 V/ y9 V- b( j) kdownFilter2 = decompose(filter2,M2);

%Second Filter
2 h/ ]: e5 x* v1 P1 ~: aupFilter3 = decompose(filter3,L3);
downFilter3 = decompose(filter3,M3);

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

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

%First Filter
( p, e: a  k: N2 o- t; q%
%%upsampling
intermediatex1 = zeros(L1,length(x));
, p; `) N% T. @/ n) [: \$ w) tx11 = zeros(L1,length(x)*L1);
" \1 J4 \6 O, o- |* ]! C+ mfor i = 11
( h% I2 {6 @% p# E/ ~/ K8 X! `intermediatex1(i, = filter(upFilter1(i,,1,x); %%do the filtering before upsampling
8 _7 u" N% |+ C& C$ X9 t' q1 xx11(i, = upsample(intermediatex1(i,:),L1); %%upsample by L, inserting L-1 zeros.
6 H4 P+ o4 f" Z3 A2 e5 Tend

X1 = zeros(1,length(x)*L1);
for i = 1:length(x)
X1(L1*(i-1)+11*(i-1)+L1) = x11(:,(i-1)*L1+1)';
# E3 t6 k+ y* L( x4 D7 T/ R2 Z3 oend
( h& e  F1 g+ x% e( O, I%%downsampling
4 c8 K/ J+ p$ x) m3 DX1_down = [X1 zeros(1,M1-mod(length(X1),M1))];
0 S* V  O7 @$ \output1 = zeros(1,length(X1_down)/M1);
0 Y( ?4 F0 K6 F5 [4 Q0 J* O- o6 [x1down = zeros(M1,length(X1_down)/M1);
  W0 r) j/ m. h- ~( B; Sx11down = x1down;
for i = 1:M1
* ?  ^! p" H' u! |& vfor k = 1:length(X1_down)/M1
5 [7 `8 \8 p- m  Ox1down(i,k) = X1_down(M1*(k-1)+i);
) @: B; q5 Q! e& j% W& V. p) T8 gend
( O$ a, u( a) P; k; `x11down(i,:) = filter(downFilter1(1,:),1,x1down(i,:));
output1 = output1 + x11down(i,:);
end

%Second Filter
%
1 o" Z8 u( w' A1 o%%upsampling
: {; p( Z, U9 p, ^5 ointermediatex2 = zeros(L2,length(output1));
9 g  `$ e6 u5 c+ C% t5 kx22 = 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.
5 `9 k) E; r+ O  T1 O! \% k, X9 Xend

X2 = zeros(1,length(output1)*L2);
8 j' l" t9 f" X, ^for i = 1:length(output1)
4 M' \5 H2 n! xX2(L2*(i-1)+1:L2*(i-1)+L2) = x22(:,(i-1)*L2+1)';
- v7 Z$ g! I6 W4 xend
# J+ v( r) P0 E6 [& @%%downsampling
X2_down = [X2 zeros(1,M2-mod(length(X2),M2))];
output2 = zeros(1,length(X2_down)/M2);
+ X) ~$ ?2 t2 H* r" cx2down = zeros(M2,length(X2_down)/M2);
x22down = x2down;
for i = 1:M2
- \4 w# q+ X% Dfor k = 1:length(X2_down)/M2
  T3 J4 o( p) X* ~x2down(i,k) = X2_down(M2*(k-1)+i);
end
0 L4 `  D) R! x8 I% C+ H# Xx22down(i,:) = filter(downFilter2(1,:),1,x2down(i,:));
output2 = output2 + x22down(i,:);
9 {5 @% L1 k/ O. Kend

%Third Filter
%
4 I& S8 c& S, @; j%%upsampling
intermediatex3 = zeros(L3,length(output2));
x33 = zeros(L3,length(output2)*L3);
, X% x, J; I: C& l+ Wfor i = 1:L3
$ Z: @0 ?) l% V9 r  lintermediatex3(i,:) = filter(upFilter3(i,:),1,output2); %%do the filtering before upsampling
3 I1 [$ |  B, @4 `2 fx33(i,:) = upsample(intermediatex3(i,:),L3); %%upsample by L, inserting L-1 zeros.
+ s( Z; f6 ~9 l0 i9 ]end

X3 = zeros(1,length(output2)*L3);
0 i. A: c; j& I- B6 ?* r( zfor i = 1:length(output2)
X3(L3*(i-1)+1:L3*(i-1)+L3) = x33(:,(i-1)*L3+1)';
& K0 p# c( t/ q) g3 tend
# ]4 j& v7 W, u1 M%%downsampling
$ w0 t1 u' j4 }8 I0 LX3_down = [X3 zeros(1,M3-mod(length(X3),M3))];
* ]! Z: [% k) r1 }" N) s5 toutput3 = zeros(1,length(X3_down)/M3);
x3down = zeros(M3,length(X3_down)/M3);
* _2 A- J8 C6 {4 ^  @x33down = x3down;
for i = 1:M3
for k = 1:length(X3_down)/M3
x3down(i,k) = X3_down(M3*(k-1)+i);
end
1 z* t0 W6 F- z# y. K/ h- K4 Wx33down(i,:) = filter(downFilter3(1,:),1,x3down(i,:));
5 k  d5 w8 j8 G9 \7 T$ f" n6 boutput3 = output3 + x33down(i,:);
end

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

nx = 1:length(output3);
2 a% `) y* @. |" W0 T' I2 hxorigin = sin(2*pi*fa*nx/96000);
error1 = xorigin - output3;
yin = [0 0 resample(x,320,147)];
% Z! o/ i; _- ?3 {  u. n- xerror2 = xorigin - yin;

figure(1)
8 r9 t# J. z5 q2 j* r0 `subplot(3,1,1);plot(xorigin);title('True Signal');
5 d8 }# @. S# u! J; Hsubplot(3,1,2);plot(output3);title('Resampled Signal');
# e. O6 a% A) w8 S  Usubplot(3,1,3);plot(error1);title('Error');

figure(2)
7 q! k2 O/ j4 J+ Q8 m2 E/ z. |subplot(3,1,1);plot(xorigin);title('True Signal');
subplot(3,1,2);plot(yin);title('Built-in Resampled Signal');
+ l' c3 q0 l, b& Qsubplot(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))];
+ S8 j+ E  z  Q& q# felse
8 P9 k% i: r; L& ?2 i- ?X = x;
end

x_decompose = zeros(factor,length(X)/factor);
for i = 1:factor
for k = 1:length(X)/factor
x_decompose(i,k) = X(factor*(k-1)+i);
end
end

end

%Filter Generation
fs = 44100;
L1 = 8;
M1 = 7;
L2 = 4;
: a# B$ ?" p1 @0 D- V5 R; t2 UM2 = 3;
L3 = 10;
" \  f1 \4 m% U# i/ P; n, EM3 = 7;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
2 m; Q/ w1 s# [; G( O%%First Filter
f1 = [20000 24000];
a1 = [1 0];
3 S: c- i% _* Bdev1 = [0.01 10^-8];
fs1 = fs*L1;
& w7 E; l( J* h/ ]% k[n1 fo1 ao1 w1] = firpmord(f1,a1,dev1,fs1);
6 t2 G+ k2 z9 ~/ r! q7 O/ U0 J6 wb1 = firpm(n1,fo1,ao1,w1);

%%Second Filter
f2 = [20160 30240];
a2 = [1 0];
dev2 = [0.01 10^-8];
fs2 = fs1*L2/M1;
3 U% k$ M; \$ Z0 c2 H[n2 fo2 ao2 w2] = firpmord(f2,a2,dev2,fs2);
6 F% {: @- r, v6 \3 Yb2 = firpm(n2,fo2,ao2,w2);

%%Third Filter
f3 = [16800 50400];
a3 = [1 0];
dev3 = [0.01 10^-8];
fs3 = fs2*L3/M2;
( k) E, s. u4 Z& i[n3 fo3 ao3 w3] = firpmord(f3,a3,dev3,fs3);
b3 = firpm(n3,fo3,ao3,w3);

木棉花_MM

这个……