Hi,
Many thanks for your suggestions.
I have replaced the abbreviations (MDCT/IMDCT) with the full names,
and added the GPL header to the codes.
Regards,
Shuhua Zhang, Aug. 1, 2009
On Sat, Aug 1, 2009 at 4:10 PM, Søren Hauberg <so...@hauberg.org> wrote:
> Hi
>
> fre, 31 07 2009 kl. 21:30 +0800, skrev Shuhua Zhang:
> > I do research in digital audio compression,
> > so times need mdct and imdct transform,
> > which is the most widely used one in audo coding,
> > such as MP3, AAC, OGG Vorbis, WMA.
> >
> > Here I post my codes for the two transforms as
> > well as a funtion to generate KBD window for AAC.
>
> This code looks good to me. I have a few comments:
>
> 1) You need to add license information to the code. If you release the
> code under the GPL that means you have to add the standard GPL header to
> each file. Similarly if you choose another license.
>
> 2) When you use abbreviations in the documentation (like MDCT) it would
> be nice if write what it stands for. That is, write
>
> % IMDCT Compute the Inverse Modified Cosine Transform (IMDCT)
> through FFT
>
> instead of
>
> % IMDCT Compute IMDCT through FFT
>
> But besides these minor issues things look good. Please fix these and
> resend the code to this list (it is really important that we make sure
> the license information is correct).
>
> Søren
>
>
function y = mdct(x)
% MDCT Compute Modified Discrete Cosine Transform (MDCT) through FFT
%
% Synopsis:
% y = mdct(x)
%
% Parameters:
% x = Input data
%
% Description:
% MDCT is the most widely used transform in audio coding, such as
% AC-3, MP3, AAC, OGG Vorbis, and WMA. For more information, see
% http://en.wikipedia.org/wiki/Modified_discrete_cosine_transform
%
% Suppose N=size(x,1)/2. The MDCT transform matrix is
% C=cos(pi/N*([0:2*N-1]'+.5+.5*N)*([0:N-1]+.5));
% then MDCT spectrum is
% y=C'*x;
%
% The fast algorithm used here consists of 4 steps:
% (1) fold column-wise 2*N rows into N rows
% (2) complex arrange the N rows into N/2 rows
% (3) pre-twiddle, N/2-CFFT, post-twiddle
% (4) reorder to form the MDCT spectrum
% In fact, (2)-(4) is a fast DCT-IV algorithm.
%
% Assumption:
% (1) x is real
% (2) size(x,1) is a multiple of 4
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%
% Author: Shuhua Zhang <shuhua.zh...@gmail.com>
% Created: 31 July 2009
nr=size(x,1);
% column oriented operation
if(nr==1)
x=x';
end
if(mod(size(x,1),4)~=0)
error(" Expecting MDCT length a multiple of 4");
end
N=size(x,1)/2;
persistent old_size w
if(isempty(old_size) || old_size~=N)
% twiddle factors
w=sparse(1:N/2, 1:N/2, exp(-j*2*pi/(2*N)*([0:N/2-1]'+.125)));
w=N^-.25*w;
old_size=N;
end
% oddly fold the first N rows
xp(1:N/2,:)=x(N/2+1:N,:)-x(N/2:-1:1,:);
% evenly fold the last N rows
xp(N/2+1:N,:)=x(N+1:3*N/2,:)+x(2*N:-1:3*N/2+1,:);
% complex arrange into N/2 rows
xc=xp(1:2:N,:)-j*xp(N:-2:1,:);
% pre-twiddle, N/2-CFFT, and post-twiddle
yc=w*fft(w*xc);
% map to MDCT
y(1:2:N,:)=imag(yc(1:N/2,:));
y(2:2:N,:)=real(yc(N/2:-1:1,:));
% row input, row output
if(nr==1)
y=y';
end
function x = imdct(y)
% IMDCT Compute Inverse Modified Discrete Cosine Transform (IMDCT) through
% IFFT
%
% Synopsis:
% x = imdct(y)
%
% Parameters:
% y = Input data
%
% Description:
% MDCT is the most widely used transform in audio coding, such as
% AC-3, MP3, AAC, OGG Vorbis, and WMA. For more information, see
% http://en.wikipedia.org/wiki/Modified_discrete_cosine_transform
%
% Suppose N=size(y,1). The MDCT transform matrix is
% C=cos(pi/N*([0:2*N-1]'+.5+.5*N)*([0:N-1]+.5));
% then given a MDCT spectrum y, the time signal x is
% x=C*y;
%
% The fast algorithm used here consists of 4 steps:
% (1) complex arrange N rows into N/2 rows
% (2) pre-twiddle, N/2-CIFFT, post-twiddle
% (3) map to N real rows
% (4) unfold the N rows into 2*N rows
% This is the inverse of the fast MDCT algorithm
%
% Assumption:
% (1) y is real
% (2) size(y,1) is a multiple of 2
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%
% Author: Shuhua Zhang <shuhua.zh...@gmail.com>
% Created: 31 July 2009
nr=size(y,1);
% column oriented operation
if(nr==1)
y=y';
end
if(mod(size(y,1),2)~=0)
error(" Expecting IMDCT length a multiple of 2");
end
N=size(y,1);
persistent old_size w
if(isempty(old_size) || old_size~=N)
% twiddle factors
w=sparse(1:N/2, 1:N/2, exp(j*2*pi/(2*N)*([0:N/2-1]'+.125)));
w=N^.25*w;
old_size=N;
end
% complex arrange into N/2 rows
yc=y(N:-2:2,:)+j*y(1:2:N,:);
% pre-twiddle, N/2-CIFFT, and post-twiddle
xc=w*ifft(w*yc);
% map to N rows
xp(1:2:N,:)=real(xc);
xp(N:-2:1,:)=-imag(xc);
% odd unfold to first N rows
x(1:N/2,:)=-xp(N/2:-1:1);
x(N/2+1:N,:)=xp(1:N/2);
% even unfold to last N rows
x(N+1:3*N/2,:)=xp(N/2+1:N,:);
x(3*N/2+1:2*N,:)=xp(N:-1:N/2+1,:);
% row input, row output
if(nr==1)
x=x';
end
function y = kbdwin(n, alpha)
% KBDWIN Generate Kaiser-Bessel-Derived (KBD) Window
%
% Synopsis:
% y = kbdwin(n, alpha)
%
% Parameters:
% n = KBD window length
% alpha = shape parameter, larger for higher sideband attenuation
%
% Description:
% Window function used in AC-3 and AAC, for more information, see
% http://en.wikipedia.org/wiki/Kaiser_window
%
% Assumption:
% (1) n is a even number
% (2) alpha is real
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%
% Author: Shuhua Zhang <shuhua.zh...@gmail.com>
% Created: 31 July 2009
% Modified Bessel function of the first kind, 0th order
x=pi*alpha*(1-(4*[0:n/2]'/n-1).^2).^.5;
x=besseli(0,x);
% accumulation
y=zeros(n,1);
y(1:n/2+1)=cumsum(x);
% normalization and mirror
y(1:n/2)=((y(1:n/2)/y(n/2+1))).^.5;
y(n/2+1:n)=y(n/2:-1:1);
------------------------------------------------------------------------------
Let Crystal Reports handle the reporting - Free Crystal Reports 2008 30-Day
trial. Simplify your report design, integration and deployment - and focus on
what you do best, core application coding. Discover what's new with
Crystal Reports now. http://p.sf.net/sfu/bobj-july
_______________________________________________
Octave-dev mailing list
Octave-dev@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/octave-dev
