Hi
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.
Shuhua Zhang,
Ph.D candidate,
Dept. of EE, Tsinghua University
function y=mdct(x)
% MDCT Compute MDCT through FFT
%
% Synopsys:
% y = mdct(x)
%
% Parameters:
% x = Input data
%
% Description:
% MDCT is the most widely used transform in audio coding, such as
% MP3, AAC, OGG Vorbis, WMA. see
% http://en.wikipedia.org/wiki/Modified_discrete_cosine_transform
% for more information.
%
% 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
%
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 IMDCT through FFT
%
% Synopsys:
% x = imdct(y)
%
% Parameters:
% y = Input data
%
% Description:
% MDCT is the most widely used transform in audio coding, such as
% MP3, AAC, OGG Vorbis, WMA. see
% http://en.wikipedia.org/wiki/Modified_discrete_cosine_transform
% for more information.
%
% 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
%
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 Window
%
% Synopsys:
% y = kbdwin(n, alpha)
%
% Parameters:
% n = KBD window length
% alpha = shape parameter, larger for higher sideband attenuation
%
% Description:
% Window function used in AAC, see
% http://en.wikipedia.org/wiki/Kaiser_window for more information
%
% Assumption:
% (1) n is a even number
% (2) alpha is real
%
% 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
