hello all this is my Java Code for Boll's 1979 Spectral Subtraction Algorithm in Java . this is based on Esfandiar Zavarehei Matlab code.
i have added some samples. the code works very good and i'm using Jtransforms lib for the FFT and IFFT
i need my code for a university project in Android. so i need to optimize it for real-time Android.
if anyone can please help. or give me a direction
thanks a lot
this is my full code http://www.4shared.com/zip/fOAyocmk/testing.html
public class GGSpecSub {
//Constructor
public GGSpecSub()
{
}
//repmat(1:W,N,1) - do this alreay transpose
public double[][] reprowtrans(int end, int replications)
{
double[][] result = new double[end +1][replications+1];
for (int x = 1; x <= end; x++) {
for (int y = 1; y <= replications; y++) {
result[x][y] = x ;
}
}
return result;
}
//repmat((0:(N-1))'*SP,1,W) this is already transposed
public double[][] repcoltrans(int end, double multiplier, int replications)
{
double[][] result = new double[replications+1][end+1];
for (int x = 1; x <= replications; x++) {
for (int y = 1; y <= end ; y++) {
result[x][y] = (y-1)*multiplier;
}
}
return result;
}
//repmat(Window,1,N)
public double[][] repvector(double[] vec, int replications)
{
double[][] result = new double[vec.length][replications];
for (int x = 0; x < vec.length; x++) {
for (int y = 0; y < replications; y++) {
result[x][y] = vec[x];
}
}
return result;
}
public double[] HammingWindow(int size)
{
double[] window = new double[size];
for (int i = 0; i < size; i++)
{
window[i] = 0.54-0.46 * (Math.cos(2.0 * Math.PI * i / (size-1)));
}
return window;
}
public double[] HanningWindow(int size)
{
double[] window = new double[size];
for (int i = 0; i < size; i++)
{
window[i] = 0.5 * (1.0 - Math.cos(2.0 * Math.PI * i / size));
}
return window;
}
public MatrixAndSegments segment (double[] signal_in,int samplesPerWindow, double shiftPercentage, double[] window)
{
//default shiftPercentage = 0.4
//default samplesPerWindow = 256 //W
//default window = hanning
int L = signal_in.length;
shiftPercentage = fix(samplesPerWindow * shiftPercentage); //SP
int numberOfSegments = fix ( (L - samplesPerWindow)/ shiftPercentage + 1); //N
double[][] reprowMatrix = reprowtrans(samplesPerWindow,numberOfSegments);
double[][] repcolMatrix = repcoltrans(numberOfSegments, shiftPercentage,samplesPerWindow );
//Index=(repmat(1:W,N,1)+repmat((0:(N-1))'*SP,1,W))';
double[][] index = new double[samplesPerWindow+1][numberOfSegments+1];
for (int x = 1; x < samplesPerWindow+1; x++ )
{
for (int y = 1 ; y < numberOfSegments + 1; y++) //numberOfSegments was 3
{
index[x][y] = reprowMatrix[x][y] + repcolMatrix[x][y];
}
}
//hamming window
double[] hammingWindow = this.HammingWindow(samplesPerWindow);
double[][] HW = repvector(hammingWindow, numberOfSegments);
double[][] seg = new double[samplesPerWindow][numberOfSegments];
for (int y = 1 ; y < numberOfSegments + 1; y++)
{
for (int x = 1; x < samplesPerWindow+1; x++)
{
seg[x-1][y-1] = signal_in[ (int)index[x][y]-1 ] * HW[x-1][y-1];
}
}
MatrixAndSegments Matrixseg = new MatrixAndSegments(numberOfSegments,seg);
return Matrixseg;
}
//not sure about casting to INT
public int fix(double val) {
if (val < 0) {
return (int) Math.ceil(val);
}
return (int) Math.floor(val);
}
public double[][] angle(double[][] array,int samplesPerWindow,int numberOfSegments)
{
//array[2*k] = Re[k],
//array[2*k+1] = Im[k]
double[][] phase = new double[fix(samplesPerWindow/2) +1][numberOfSegments];
if (fix(samplesPerWindow*2/2) + 2 % 2 != 0)
{
for (int y = 0; y < numberOfSegments; y++ )
{
for (int x = 0, k = 0; x <= fix(samplesPerWindow*2/2); x = x + 2, k++)
{
phase[k][y] = Math.atan2(array[x+1][y],array[x][y]);
}
}
}
else
{
for (int y = 0; y < numberOfSegments; y++ )
{
for (int x = 0, k = 0; x < fix(samplesPerWindow*2/2) + 2; x = x + 2, k++)
{
phase[k][y] = Math.atan2(array[x+1][y],array[x][y]);
}
}
}
return phase;
}
public double[][] complexAbs(double[][] array,int samplesPerWindow,int numberOfSegments)
{
//array[2*k] = Re[k],
//array[2*k+1] = Im[k]
double[][] abs = new double[fix(samplesPerWindow/2) +1][numberOfSegments];
if (fix(samplesPerWindow*2/2) + 2 % 2 != 0)
{
for (int y = 0; y < numberOfSegments; y++ )
{
for (int x = 0, k = 0; x <= fix(samplesPerWindow*2/2); x = x + 2, k++)
{
abs[k][y] = Math.sqrt( ((array[x+1][y])*(array[x+1][y])) + ((array[x][y])*(array[x][y])) );
}
}
}
else
{
for (int y = 0; y < numberOfSegments; y++ )
{
for (int x = 0, k = 0; x < fix(samplesPerWindow*2/2); x = x + 2, k++)
{
abs[k][y] = Math.sqrt( ((array[x+1][y])*(array[x+1][y])) + ((array[x][y])*(array[x][y])) );
}
}
}
return abs;
}
//does mean of each row and divide by the number of cols
public double[] mean(double[][] array,int samplesPerWindow,int numberOfSegments,int nis )
{
double[] meanArray = new double[samplesPerWindow];
for (int x = 0 ; x < samplesPerWindow ; x++)
{
for (int y = 0; y < nis ; y++)
{
meanArray[x] += array[x][y];
}
meanArray[x] = meanArray[x] / nis;
}
return meanArray;
}
public vadResult vad(double[] signal,double[] noise,int noiseCounter)
{
vadResult vadRes = new vadResult();
int noiseMargin = 3;
int hangover = 8;
double sumSpectralDist = 0.0;
int FreqResol = signal.length; //FreqResol=length(signal);
double[] spectralDist = new double[FreqResol];
//SpectralDist= 20*(log10(signal)-log10(noise));
for(int i = 0; i < signal.length ; i++){
double temp = 20 * (Math.log10(signal[i]) - Math.log10(noise[i]));
spectralDist[i] = (temp < 0) ? 0 : temp;
sumSpectralDist += spectralDist[i];
}
double dist = sumSpectralDist/(double)signal.length;
if (dist < noiseMargin){
vadRes.noiseFlag = true;
noiseCounter = noiseCounter + 1;
}
else{
vadRes.noiseFlag = false;
noiseCounter = 0;
}
if (noiseCounter > hangover){
vadRes.speechFlag = false;
}
else{
vadRes.speechFlag = true;
}
vadRes.noiseCounter = noiseCounter;
return vadRes;
}
public double[] getColumn(double[][] signal,int i)
{
double[] column = new double[signal.length];
for (int j = 0 ; j < signal.length ; j++)
{
column[j] = signal[j][i];
}
return column;
}
public double[] OverLapAdd2(double[][] XNEW,int XNEWcolNum, double[][]yphase,int windowLen, int shiftLen)
{
//double[] res= new double[XNEWcolNum];//check me
double[][] Spec = new double[2*XNEW.length][XNEWcolNum];// x2 fro complex numbers
if (fix(shiftLen) != shiftLen)
shiftLen = fix(shiftLen);
int FreqRes = XNEW.length;//number of Rows in Xnew
int FrameNum = XNEWcolNum; // number of cols in XNEW
//array[2*k] = Re[k],
//array[2*k+1] = Im[k]
for (int y = 0; y < XNEWcolNum; y++ )
{
for (int x = 0,k = 0 ; x < 2*XNEW.length; x=x+2, k++)
{
double length = Math.exp(0);
Spec[x][y] =XNEW[k][y]*length*Math.cos(yphase[k][y]); //real
Spec[x + 1][y] =XNEW[k][y]*length*Math.sin(yphase[k][y]); //im
//Spec=XNEW.*exp(j*yphase);
}
}
double[][] flipudMatrix= this.flipud(Spec,XNEWcolNum,windowLen);//if mod(windowLen,2) %if FreqResol is odd
// Spec=[Spec;flipud(conj(Spec(2:end,:)))];
//else
// Spec=[Spec;flipud(conj(Spec(2:end-1,:)))];
//end
//we do NOT join the Spec and FlipUdmatrix it's a waste of time
double[] signal_out = new double[ (((FrameNum -1) * (shiftLen)) + windowLen)];
DoubleFFT_1D fft = new DoubleFFT_1D((flipudMatrix.length + Spec.length )/2);
for (int i = 0 ; i < FrameNum; i++)
{
int start = ((i /*- 1*/) * shiftLen) + 1; //start=(i-1)*ShiftLen+1;
double[] specAndFlipudcol = new double[flipudMatrix.length + Spec.length];
for (int j = 0 ; j < (Spec.length) ; j++)
{
specAndFlipudcol[j] = Spec[j][i];
}
for (int j = 0; j < flipudMatrix.length ;j++)
{
specAndFlipudcol[j + Spec.length] = flipudMatrix[j][i]; // spec=Spec(:,i);
}
fft.realInverse(specAndFlipudcol,true);
for (int j = start,k = 0 ; j < start + windowLen - 1 ;k++, j++)
signal_out[j] = signal_out[j] + specAndFlipudcol[k];
int gili = 0;
}
return signal_out;
//return res;
}
public double[] SSBoll79(double[] signal_in,int fs)
{
int length = signal_in.length;
int numberOfSegments;
int samplesPerWindow = fix(0.025*fs); //W
double shiftPercentage = 0.4;
double IS=0.25;
int nis = (int)this.fix((IS*fs-samplesPerWindow)/(shiftPercentage * samplesPerWindow)+1); //%number of initial silence segments
double[] window = new double[samplesPerWindow]; //wnd = hamming(W)
//need to return both the Matrix and number of segments
MatrixAndSegments Matrixres = this.segment(signal_in, samplesPerWindow, shiftPercentage, window); //y=segment(signal,W,SP,wnd);
double [][] newRes = new double[samplesPerWindow*2][Matrixres.numberOfSegments];
double [] colForFFT = new double [samplesPerWindow*2];
DoubleFFT_1D fft = new DoubleFFT_1D(samplesPerWindow);
for(int y = 0; y < Matrixres.numberOfSegments; y++)
{
//copy the original col into a col and and a col of zeros before FFT
for(int x = 0; x < samplesPerWindow; x++)
{
colForFFT[x] = Matrixres.res[x][y];
}
//fft on each col of the matrix
fft.realForwardFull(colForFFT); //Y=fft(y,nfft);
//copy the output of col*2 size into a new matrix
for(int x = 0; x < samplesPerWindow*2; x++)
{
newRes[x][y] = colForFFT[x];
}
}
double[][] yphase = this.angle(newRes,samplesPerWindow,Matrixres.numberOfSegments); //YPhase=angle(Y(1:fix(end/2)+1,:)); %Noisy Speech Phase
newRes = this.complexAbs(newRes, samplesPerWindow, Matrixres.numberOfSegments); //Y=abs(Y(1:fix(end/2)+1,:)).^Gamma;%Specrogram
int numberOfFrames = Matrixres.numberOfSegments; //numberOfFrames=size(Y,2); //number of cols in newRes matrix
int FreqResol = (this.fix(samplesPerWindow) + 2) / 2; //FreqResol=size(Y,1); //number of rows in newRes matrix
double[] N = this.mean(newRes, FreqResol, Matrixres.numberOfSegments,nis); //N=mean(Y(:,1:NIS)')'; %initial Noise Power Spectrum mean
double[] NRM = new double[N.length]; // NRM=zeros(size(N));% Noise Residual Maximum (Initialization)
int noiseCounter = 0 ;
int noiseLength = 9;
double beta = 0.03;
double[][] YS = new double[newRes.length][];
for(int i = 0; i < newRes.length; i++)
YS[i] = newRes[i].clone(); //YS=Y; %Y Magnitude Averaged
/*MATLAB CODE
* for i=2:(numberOfFrames-1)
YS(:,i)=(Y(:,i-1)+Y(:,i)+Y(:,i+1))/3;
end
*/
for (int y = 1; y < numberOfFrames-1; y++)
{
for (int x = 0; x < YS.length; x++ )
{
if (y==0){}
else {
YS[x][y] = (newRes[x][y-1] + newRes[x][y] + newRes[x][y+1]) / 3.0; //can be optimized Gilad
}
}
}
double[][] X = new double[NRM.length][numberOfFrames];
for (int i = 0; i < numberOfFrames; i++)
{
double[] YScolumn = this.getColumn(YS, i);
double[] column = this.getColumn(newRes, i);
vadResult tempVadResult = this.vad(column, N, noiseCounter); //[NoiseFlag, SpeechFlag, NoiseCounter, Dist]=vad(Y(:,i).^(1/Gamma),N.^(1/Gamma),NoiseCounter); %Magnitude Spectrum Distance VAD
noiseCounter = tempVadResult.noiseCounter; // keep the result of the noise counter from last loop
if (tempVadResult.speechFlag == false)
{ //if SpeechFlag==0
//N=(NoiseLength*N+Y(:,i))/(NoiseLength+1); %Update and smooth noise
for (int j = 0; j < N.length ; j++)
{
N[j] = (noiseLength * N[j] + column[j])/ (noiseLength + 1);
YScolumn[j] = YScolumn[j]- N[j];
NRM[j] = Math.max(NRM[j], YScolumn[j]); //NRM=max(NRM,YS(:,i)-N);%Update Maximum Noise Residue
X[j][i] = beta * column[j]; // X(:,i)=Beta*Y(:,i);
}
}
//speechFlag == 1
else
{
double[] D = new double[YScolumn.length];
for (int j = 0 ;j < YScolumn.length ; j++)
{
D[j] = YScolumn[j]- N[j];
}
if ((i > 0) && (i < numberOfFrames-1))
{ /* for j=1:length(D)
if D(j)<NRM(j)
D(j)=min([D(j) YS(j,i-1)-N(j) YS(j,i+1)-N(j)]);
end
end
*/
for (int j = 0 ; j < D.length ; j++)
{
if (D[j] < NRM[j])
{
D[j] = Math.min( Math.min(D[j], YS[j][i-1] - N[j]), (YS[j][i+1]-N[j]) );
}
}
}
for (int j = 0; j <D.length;j++)
{
X[j][i] = (D[j] < 0.0) ? 0: D[j]; //X(:,i)=max(D,0);
}
}
}
double[] output= this.OverLapAdd2(X, /*XNEWcolNum*/ numberOfFrames, yphase, /*windowLen*/samplesPerWindow, (int)(shiftPercentage*samplesPerWindow));
return output;
}
public double[][] flipud(double[][] matrix,int columns,int windowLen)
{
double[][] conjMatrix = new double[matrix.length][columns]; // does conj to all the complex numbers
int end;
conjMatrix = conj(matrix,columns);
double[][] flipudMatrix;
if ((windowLen) % 2 == 0){
end = matrix.length - 4;
flipudMatrix = new double[end][columns];
}
else{
end = matrix.length - 2;
flipudMatrix = new double[end][columns];
}
for (int i = end,k = 0 ; i >= 2 ;i = i-2, k= k + 2)
{
for(int j = 0 ; j < columns ; j++)
{
flipudMatrix[k][j] = conjMatrix[i][j];//real part
flipudMatrix[k+1][j] = conjMatrix[i+1][j];//imaginary part
}
}
return flipudMatrix;
}
public double[][] conj(double[][] originalMatrix,int columns)
{
double [][] matrix = new double[originalMatrix.length][];
for(int i = 0; i < originalMatrix.length; i++)
matrix[i] = originalMatrix[i].clone();
int end = ((matrix.length/2) % 2 == 0)? matrix.length-2 :matrix.length ;
for (int i = 3; i < end ; i=i+2)//we need just the row that are odd and we don't change the first 2 lines // so no 0+1 lines , line 2 is real . line 3 is the first img we need
{
for (int j = 0 ; j < columns ; j++)
{
matrix[i][j] = -1*matrix[i][j];
}
}
return matrix;
}
}
