[Scilab-users] [Users] Discrete cosine transform

[Serge Steer]

With Scilab 5.4.0 Windows versions using the MKL library, the dct and 
dst function does not work (the associated routines are not present in 
the MKL library)
One solution is to use a non MKL version
You can also use the attached Scilab function based on fft (slower)

Serge Steer
Le 26/03/2013 09:36, TViT a écrit :
> Say please why the example, given in the help does not work? What
> communication(connection) DCT and FFTw?
> In my representation FFT through one script is carried out, and DCT through
> another script however dct.sci has not found in Scilab 5.4.0
>
> //----------------------------------
> -->d=dct(s);
>           !--error 999
> dct: Creation of requested fftw plan failed.
>
> //----------------------------------
> // build a sampled at 1000hz  containing pure frequencies
> // at 50 and 70 Hz
> sample_rate=1000;
> t = 0:1/sample_rate:0.6;
> N=size(t,'*'); //number of samples
> s=sin(2*%pi*50*t)+sin(2*%pi*70*t+%pi/4)+grand(1,N,'nor',0,1);
> d=dct(s);
> // zero low energy components
> d(abs(d)<1)=0;
> size(find(y1<>0),'*') //only 30 non zero coefficients out of 600
> clf;plot(s,'b'),plot(dct(d,1),'r')
>
>
>
> --
> View this message in context: http://mailinglists.scilab.org/Users-Discrete-cosine-transform-tp4024712p4026359.html
> Sent from the Scilab users - Mailing Lists Archives mailing list archive at Nabble.com.
> _______________________________________________
> users mailing list
> users at lists.scilab.org
> http://lists.scilab.org/mailman/listinfo/users
>

-------------- next part --------------
function y=dct(A,varargin)
  if A==[] then y=A;return;end
  if size(A,'*')==1 then y=A;return;end
  if size(varargin)==0 then
    d=-1
  else
    d=varargin(1)
    varargin(1)=null()
  end
  
  select size(varargin)
  case 0 then //dct(A [,sign])
     y=dct_nd(A,d)
  case 1 then //dct(A [,sign],selection)
  case 3 then //dct(A [,sign],dims,incrs)
  end
endfunction

function y= dct_nd(y,flag)
  for k=1:ndims(y)
    y=dct_1d(y,flag,k)
  end
endfunction

function y = dct_1d (x,flag,sel)
//Ref:  A. K. Jain, "Fundamentals of Digital Image Processing", pp. 150-153.
  if x==[] then y=x,return,end
  if argn(2)==1 then flag=-1,end
  szx=size(x);
  if argn(2)<3 then 
    if size(find(szx>1),'*')==1 then sel=k;else sel=1;end
  end
  n=szx(sel)
  if n==1 then y=x,return,end
  if sel<>1 then
    sz1=1:ndims(x);sz1([1 sel])=sz1([sel 1]);
    x=permute(x,sz1) 
    szx=permute(szx,sz1)
  end
  x=matrix(x,n,-1);
  m=size(x,2);
  if flag==-1 then
    w=[sqrt(0.25/n);sqrt(0.5/n)*exp((-%i*%pi/2/n)*(1:n-1)')]*ones(1,m);
    if isreal(x)&modulo(n,2)==0 then
      y=2*real(w.*fft([x(1:2:n,:); x(n:-2:2,:)],flag,1));
    else
      y=fft([x;x($:-1:1,:)],flag,1);
      y=w.*y(1:n,:);
      if isreal(x) then y=real(y); end
    end
  else //inverse transform
    if isreal(x)&modulo(n,2)==0 then
      w = [sqrt(n/4);sqrt(n/2)*exp((%i*%pi/2/n)*(1:n-1)')]*ones(1,m);
      y=fft(w.*x,flag,1);
      y([1:2:n,n:-2:1],:)=2*real(y);
    else
      w = [sqrt(4*n);sqrt(2*n)*exp((%i*%pi/2/n)*(1:n-1)')]*ones(1,m);
      y=[x.*w;zeros(1,nc);-%i*w(2:$,:).*x($:-1:2,:)];
      y=fft([x.*w;zeros(1,nc);-%i*w(2:$,:).*x($:-1:2,:)],flag,1);
      y=y(1:n,:);
      if isreal(x) then y=real(y); end
    end
  end
  if sel<>1 then  y=permute(y,sz1) ;end
  y=matrix(y,szx)

endfunction

function y=idct(x)
  y=dct(x,1)
endfunction
function y=dst(x,flag)
  [mx,nx]=size(x);
  if argn(2)==1 then flag=-1,end
  if size(x,1)==1 then 
    //transform x as a column vector
    x=matrix(x,-1,1);
    n=mx*nx;
    nc=1;
  else
    n=mx;
    nc=nx
  end
  if n==1 then y=matrix(x,mx,nx),return,end
  y=fft([zeros(1,nc);x;zeros(1,nc);-x($:-1:1,:)],-1,1)
  if isreal(x) then
    y=-imag(y(2:n+1,:))/2
  else
    y=y(2:n+1,:)/(-2*%i)
  end
  if flag==1 then y=2/(n+1)*y;end
  y=matrix(y,mx,nx);
endfunction
function y=idst(x)
  y=dst(x,1)
endfunction