|
0
|
1 ## Copyright (C) 2000 Paul Kienzle |
|
|
2 ## |
|
|
3 ## This program is free software; you can redistribute it and/or modify |
|
|
4 ## it under the terms of the GNU General Public License as published by |
|
|
5 ## the Free Software Foundation; either version 2 of the License, or |
|
|
6 ## (at your option) any later version. |
|
|
7 ## |
|
|
8 ## This program is distributed in the hope that it will be useful, |
|
|
9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
|
|
11 ## GNU General Public License for more details. |
|
|
12 ## |
|
|
13 ## You should have received a copy of the GNU General Public License |
|
|
14 ## along with this program; if not, write to the Free Software |
|
|
15 ## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA |
|
|
16 |
|
|
17 ## usage: y=resample(x,p,q,d) |
|
|
18 ## |
|
|
19 ## Change the sample rate of x by a factor of p/q. Note that p and q do |
|
|
20 ## not need to be integers since this routine does not use a polyphase |
|
|
21 ## rate change algorithm, but instead uses bandlimited interpolation, |
|
|
22 ## wherein the continous time signal is estimated by summing the sinc |
|
|
23 ## functions of the nearest neighbouring points up to distance d. |
|
|
24 ## |
|
|
25 ## This is discussed in: |
|
|
26 ## J. O. Smith and P. Gossett (1984). A flexible sampling-rate |
|
|
27 ## conversion method. In ICASSP-84, Volume II, pp. 19.4.1-19.4.2. |
|
|
28 ## New York: IEEE Press. |
|
|
29 ## See the authors page at: http://www-ccrma.stanford.edu/~jos/resample/ |
|
|
30 ## |
|
|
31 ## Note that the resampling is not yet very fast or very good, but it is |
|
|
32 ## very flexible. |
|
|
33 ## |
|
|
34 ## Example |
|
|
35 ## ## Speech example |
|
|
36 ## [x, fs] = auload(file_in_loadpath("sample.wav")); |
|
|
37 ## sound(resample(x,16000,fs), 16000); # resample at 16 kHz |
|
|
38 ## |
|
|
39 ## ## Example from interp1 |
|
|
40 ## xf=0:0.05:10.95; yf = sin(2*pi*xf/5); |
|
|
41 ## xp=0:10; yp = sin(2*pi*xp/5); |
|
|
42 ## r = resample(yp,xp(2),xf(2)); |
|
|
43 ## plot(xf,yf,';original;',xf,r,';resample;',xp,yp,'*;;'); |
|
|
44 ## |
|
|
45 ## Note that resample computes all samples up to but not including time |
|
|
46 ## n+1. If you are increasing the sample rate, this means that it will |
|
|
47 ## generate samples beyond the end of the time range of the original |
|
|
48 ## signal. That is why xf must goes all the way to 10.95 in the example. |
|
|
49 |
|
|
50 ## TODO: Fix so that audible clicking goes away. |
|
|
51 ## TODO: Change to a faster algorithm. |
|
|
52 ## TODO: Test on a chirp signal. |
|
|
53 |
|
|
54 function y=resample(x,p,q,order,beta) |
|
|
55 if (nargin < 2 || nargin > 5) |
|
|
56 usage("y=resample(x,p,q,order)"); |
|
|
57 endif |
|
|
58 |
|
|
59 if (nargin < 3), q=1; endif |
|
|
60 if (nargin < 4), order = 5; endif |
|
|
61 |
|
|
62 ## ## chain to decimate/interpolate if appropriate |
|
|
63 ## if p==1 && q==fix(q) |
|
|
64 ## y=decimate(x,q); order? |
|
|
65 ## return; |
|
|
66 ## elseif q==1 && p==fix(p) |
|
|
67 ## y=interp(x,q); order? |
|
|
68 ## return; |
|
|
69 ## endif |
|
|
70 |
|
|
71 transpose = rows(x)==1; |
|
|
72 if transpose, x = x.'; endif |
|
|
73 |
|
|
74 ## if rate reduction, apply antialiasing filter first |
|
|
75 r=p/q; |
|
|
76 if (r < 1) |
|
|
77 b = fir1(2*order+1, r); |
|
|
78 x = fftfilt(b, x); |
|
|
79 endif |
|
|
80 |
|
|
81 ## Determine the new sampling times, and their distance to the old |
|
|
82 ## ones. Note that the new series should be the maximum that can |
|
|
83 ## be contained in the old series without going over the time |
|
|
84 ## allotted to the old series. In short, you have to go a little |
|
|
85 ## beyond the last sample of the old series if your new sampling |
|
|
86 ## rate is higher. |
|
|
87 t=[1:1/r:length(x)+1-1/r]'; # the sampling points of the new series |
|
|
88 idx = fix(t); # the nearest old point |
|
|
89 t = t-idx; # distance to the nearest old point |
|
|
90 |
|
|
91 ## generate the new series by summing the sinc functions of the |
|
|
92 ## nearest neighbour points implicit in the continuous time |
|
|
93 ## expansion of the old series. This new series is truncated |
|
|
94 ## to +/- order nearest neighbours. For convenience, the original |
|
|
95 ## series is zero-padded before and after, implicitly setting the |
|
|
96 ## neighbours at the start of the signal to zero. |
|
|
97 x = [zeros(order,columns(x)) ; x ; zeros(order,columns(x))]; |
|
|
98 y = zeros(length(idx),columns(x)); # the new series |
|
|
99 for i=-order:order |
|
|
100 w = sinc(t-i).*(0.5+0.5*cos(pi*(t-i)/(order+0.5))); # hanning window |
|
|
101 y=y + x(idx+i+order,:).*w(:,ones(size(x,2),1)); |
|
|
102 endfor |
|
|
103 |
|
|
104 if transpose, y=y.'; endif |
|
|
105 endfunction |
|
|
106 |
|
|
107 %!demo |
|
|
108 %! xf=0:0.05:10.95; yf = sin(2*pi*xf/5); |
|
|
109 %! xp=0:10; yp = sin(2*pi*xp/5); |
|
|
110 %! r = resample(yp,xp(2),xf(2)); |
|
|
111 %! oneplot(); |
|
|
112 %! title("confirm that the resampled function matches the original"); |
|
|
113 %! plot(xf,yf,';original;',... |
|
|
114 %! xf,r(1:length(xf)),';resample;',... |
|
|
115 %! xp,yp,'*;;'); |
|
|
116 %! title(""); |
|
|
117 %! [x, fs] = auload(file_in_loadpath("sample.wav")); |
|
|
118 %! sound(resample(x,16000,fs), 16000); |
