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

loveeatmore

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

%input data
fa = 8000; %%signal frequency
fs = 44100; %%44.1kHz sampling frequency
' U1 y0 j+ ~  d$ k" Cn = 1:64;
; z2 Y% s$ Q0 S8 s7 H; M4 @x = sin(2*pi*fa*n/fs);
lengthx = length(x);
' p! K  I. B& `8 y7 B* P. \6 N! dt = n*1000/128/fs;
% A7 _# ~. @( S+ m, q; |%plot(t,x);xlabel('time in ms');

% [y fs nbits] = wavread('acoustic.wav');
2 a  X" u% h/ k  N, V  U) T% x = y';
5 q. D1 ?2 C% j; V/ {6 o& G+ E+ ~% lx = length(x);
% J8 E% s3 k0 `%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
7 f$ K% f/ s5 T- @%Resample parameters
# N( T) b" |2 Q, w- s9 k5 wL1 = 8;
* ~/ @! R+ I6 o6 qM1 = 7;
9 c- p& t9 r$ M( _. NL2 = 4;
M2 = 3;
L3 = 10;
M3 = 7;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! t5 [( i/ V2 }, v%%filters
, R3 c7 u# ^2 W9 s. r0 Iload filter.mat
filter1 = b1;
filter2 = b2;
/ w7 f% i# R+ B8 p" Q( hfilter3 = b3;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  T& y1 W' x! n3 B  n2 D% H3 W- P%Initialize the subfilters

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

%Second Filter
8 a  g) u7 ^) M" P$ YupFilter2 = decompose(filter2,L2);
downFilter2 = decompose(filter2,M2);

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

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

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

%First Filter
5 t$ Z4 L2 e4 D$ U6 P%
%%upsampling
! p0 E7 }  A: A( L, N) Aintermediatex1 = zeros(L1,length(x));
  b& o3 m/ C/ T5 px11 = zeros(L1,length(x)*L1);
for 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.
end

X1 = zeros(1,length(x)*L1);
; ~& e3 w* z, z5 G# d/ Jfor i = 1:length(x)
- \/ K0 n! j& r/ kX1(L1*(i-1)+11*(i-1)+L1) = x11(:,(i-1)*L1+1)';
" J5 ~8 q$ s# ^4 I0 B+ Aend
* h% p" U) \& k* T%%downsampling
X1_down = [X1 zeros(1,M1-mod(length(X1),M1))];
7 k3 `1 w' u/ V; f4 c! n" Poutput1 = zeros(1,length(X1_down)/M1);
+ }* K/ M  h- Lx1down = zeros(M1,length(X1_down)/M1);
x11down = x1down;
& A: T- [4 Z% u9 k% b/ y% ~for i = 1:M1
4 x- l" u: R0 n! g1 Z% \, [for k = 1:length(X1_down)/M1
x1down(i,k) = X1_down(M1*(k-1)+i);
0 L& [$ A. i; O2 Cend
x11down(i,:) = filter(downFilter1(1,:),1,x1down(i,:));
" r3 Q3 [. E- foutput1 = output1 + x11down(i,:);
end

%Second Filter
%
" a- ?! X9 x% L%%upsampling
6 e, E* _( d1 C2 S  t9 F& ^& [intermediatex2 = zeros(L2,length(output1));
x22 = zeros(L2,length(output1)*L2);
for i = 12
4 @6 D/ G6 t6 u) J: Lintermediatex2(i,:) = filter(upFilter2(i,:),1,output1); %%do the filtering before upsampling
x22(i,:) = upsample(intermediatex2(i,:),L2); %%upsample by L, inserting L-1 zeros.
( e& k; ~. _+ a& Mend

X2 = zeros(1,length(output1)*L2);
+ k) D4 b7 e4 y7 w* ]% Mfor i = 1:length(output1)
X2(L2*(i-1)+1:L2*(i-1)+L2) = x22(:,(i-1)*L2+1)';
end
%%downsampling
4 Y7 u7 E+ v; N4 D8 o, ?7 hX2_down = [X2 zeros(1,M2-mod(length(X2),M2))];
output2 = zeros(1,length(X2_down)/M2);
x2down = zeros(M2,length(X2_down)/M2);
x22down = x2down;
! r( d, h0 G1 z( t1 mfor i = 1:M2
for k = 1:length(X2_down)/M2
7 y7 C# `( k( _' |2 U' Fx2down(i,k) = X2_down(M2*(k-1)+i);
! |$ B7 Z* v: cend
. n* r" F2 P  L7 x; [! |: J# v3 Vx22down(i,:) = filter(downFilter2(1,:),1,x2down(i,:));
. |) ~) k" y* c  o: @0 L" ?/ {$ aoutput2 = output2 + x22down(i,:);
end

%Third Filter
%
: W( l. R5 [7 `+ w: {# z3 T+ r%%upsampling
) g) t/ ]9 m* @; W% I' [2 pintermediatex3 = zeros(L3,length(output2));
x33 = zeros(L3,length(output2)*L3);
' B$ j* j! j; A9 q+ S" z6 Ifor i = 1:L3
# O7 R3 G5 H, Y6 x* X" a8 ?intermediatex3(i,:) = filter(upFilter3(i,:),1,output2); %%do the filtering before upsampling
4 Q# I# ]+ Q4 Z7 G. Px33(i,:) = upsample(intermediatex3(i,:),L3); %%upsample by L, inserting L-1 zeros.
7 O) ^& r6 G' @) V5 `! {end

X3 = zeros(1,length(output2)*L3);
* J4 o* Q( }& c7 _for i = 1:length(output2)
) ~8 I) Z$ c& g8 nX3(L3*(i-1)+1:L3*(i-1)+L3) = x33(:,(i-1)*L3+1)';
end
%%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;
for i = 1:M3
for k = 1:length(X3_down)/M3
( t. T2 s- Z/ J: t8 _x3down(i,k) = X3_down(M3*(k-1)+i);
end
x33down(i,:) = filter(downFilter3(1,:),1,x3down(i,:));
output3 = output3 + x33down(i,:);
end

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

nx = 1:length(output3);
$ @  B* z: Z9 l- r& }/ V2 M; \" wxorigin = sin(2*pi*fa*nx/96000);
2 P7 O  f. c( p0 x1 h' M: b& n- eerror1 = xorigin - output3;
5 g& O6 ~6 C# k# L- eyin = [0 0 resample(x,320,147)];
/ t  e& m- {4 zerror2 = xorigin - yin;

figure(1)
subplot(3,1,1);plot(xorigin);title('True Signal');
- ]1 f1 @4 P& }4 _- y& b" |* v2 j; r2 Ksubplot(3,1,2);plot(output3);title('Resampled Signal');
subplot(3,1,3);plot(error1);title('Error');

figure(2)
$ |5 S7 @( l4 `* f- Q. ~1 x  }subplot(3,1,1);plot(xorigin);title('True Signal');
0 z8 w( f) X/ Y' N: c9 ysubplot(3,1,2);plot(yin);title('Built-in Resampled Signal');
! K7 r; T: U8 c8 Lsubplot(3,1,3);plot(error2);title('Error');

---

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

if mod(lx,factor) ~= 0;
! }( ]9 T& i# O9 uX = [x zeros(1,factor-mod(lx,factor))];
3 d. W6 A6 {: l$ U+ j9 N: Gelse
) Z. p9 _  H6 d; U8 pX = x;
end

x_decompose = zeros(factor,length(X)/factor);
2 `( T2 n$ Y3 J! ifor i = 1:factor
for k = 1:length(X)/factor
* t7 P- p/ e' B. {5 t4 `7 `( |x_decompose(i,k) = X(factor*(k-1)+i);
5 M% q6 r2 }0 K3 bend
+ F2 [  W& q7 e8 d0 gend

end

%Filter Generation
; ?2 m' f% S4 cfs = 44100;
L1 = 8;
M1 = 7;
7 B* K. i# O  `' }7 Q0 r* q/ N8 Z5 N, e8 BL2 = 4;
M2 = 3;
L3 = 10;
% L4 U( ~. L; ?. e. H  U. zM3 = 7;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- Z$ u" y" @- t$ O; l%
% U9 {# S5 {) A, o%%First Filter
: T  ?" n* F& a4 o. [1 i2 Nf1 = [20000 24000];
a1 = [1 0];
1 W4 G, |" K; C8 `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
/ u: F% V) R3 t! }# A7 Pf2 = [20160 30240];
a2 = [1 0];
dev2 = [0.01 10^-8];
+ z6 g' {6 l7 F2 P1 `. H4 D6 ]8 N" Q, Mfs2 = 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];
7 A4 _5 {* e' `dev3 = [0.01 10^-8];
3 R0 F! i; y5 |" qfs3 = fs2*L3/M2;
[n3 fo3 ao3 w3] = firpmord(f3,a3,dev3,fs3);
, S# `% L9 \  K$ X% [. Qb3 = firpm(n3,fo3,ao3,w3);

! F( [6 P; n: ]$ B

木棉花_MM

这个……