indrajith/package-zearch-temp
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
master
package-zearch-temp/UOS/examples/uos_dsp_noiseremoval.pas
Indrajith K L 6a1d80d3b8 Implements Module Player working in Form 
* Implements Embedding .xm music in Final Application
2025-06-02 15:15:16 +05:30

1203 lines
31 KiB
ObjectPascal
Raw Permalink Blame History

{This unit is part of United Openlibraries of Sound (uos)}
{
audacity which this file was converted from is GPLv2
http://www.audacityteam.org/about/license/
Converted By Andrew Haines
License : modified LGPL.
Fred van Stappen fiens@hotmail.com
}
unit uos_dsp_noiseremoval;
{$mode objfpc}{$H+}
interface
uses
Classes, SysUtils, uos_dsp_utils ;
type
PFFT = ^TFFT;
{ TFFT }
TFFT = object
BitReversed: PInteger;
SinTable: PSingle;
Points: Integer;
FPCSinTable: array of Single;
function InitializeFFT(FFTLen: Integer): PFFT; static;
procedure EndFFT;
function GetFFT(FFTLen: Integer): PFFT; static;
procedure ReleaseFFT;
procedure InverseRealFFTf(buffer: PSingle);
procedure CleanupFFT; static; // ???
procedure RealFFTf(buffer: PSingle);
procedure ReorderToTime(Buffer: PSingle; TimeOut: PSingle);
procedure ReorderToFreq(Buffer: PSingle; RealOut: PSingle; ImagOut: PSingle);
end;
type
PPSingle=^PSingle;
TSingleArray = array of Single;
TNoiseRemoval = class;
TNoiseWriteProc = procedure(ASender: TObject; AData: PSingle; ASampleCount: Integer) of Object;
{ TNoiseRemoval }
TNoiseRemoval = class
private
FDoProfile: Boolean;
FHasProfile: Boolean;
// Parameters chosen before the first phase
FSampleRate: Double;
FWindowSize: Integer;
FSpectrumSize: Integer;
FMinSignalTime: Single; // in secs
// The frequency-indexed noise threshold derived during the first
// phase of analysis
FNoiseThreshold: array of Single; // length in FSpectrumSize
// Parameters that affect the noise removal, regardless of how the
// noise profile was extracted
FSensitivity: Double;
FFreqSmoothingHz: Double;
FNoiseGain: Double; // in dB, should be negative
FAttackDecayTime: Double;
FbLeaveNoise: Boolean;
// change this later
procedure Initialize;
procedure Reset; // StartNewTrack
procedure ProcessSamples(len: Integer; Buffer:PSingle);
procedure FillFirstHistoryWindow;
procedure ApplyFreqSmoothing(ASpec: PSingle);
procedure GetProfile;
procedure RemoveNoise;
procedure RotateHistoryWindows;
procedure FinishTrack;
procedure Cleanup;
private
//FOutputTrack: PSingle; // WaveTrack;
FInSampleCount: Integer;
FOutSampleCount: Integer;
FInputPos: Integer;
FFFT: PFFT;
FFFTBuffer: PSingle; // FWindowSize
FWindow: PSingle; // FWindowSize
FFreqSmoothingBins: Integer;
FAttackDecayBlocks: Integer;
FOneBlockAttackDecay: Single;
FNoiseAttenFactor: Single;
FSensitivityFactor: Single;
FMinSignalBlocks: Integer;
FHistoryLen: Integer;
FInWaveBuffer: PSingle; // FWindowSize
FOutOverlapBuffer: PSingle; // FWindowSize
FSpectrums: array of PSingle; // FHistoryLen x FSpectrumSize
FGains: array of PSingle; // FHistoryLen x FSpectrumSize
FRealFFTs: array of PSingle; // FHistoryLen x FWindowSize
FImagFFTs: array of PSingle; // FHistoryLen x FWindowSize
FWriteProc: TNoiseWriteProc;
FInited: Boolean;
FTotalRead: QWord;
function GetNoiseProfile: TSingleArray;
procedure SetAttackDecayTime(AValue: Double);
procedure SetFreqSmoothingHz(AValue: Double);
procedure SetGain(AValue: Double);
procedure SetNoiseProfile(AValue: TSingleArray);
procedure SetSensitivity(AValue: Double);
public
constructor Create;
destructor Destroy; override;
function Init(ASampleRate: Integer): Boolean;
function Process(AData: PSingle; ASampleCount: Integer; AGetNoiseProfile: Boolean; AMoreComing: Boolean = False): Boolean;
procedure Flush; // finish writing out data in buffers.
property NoiseProfile: TSingleArray read GetNoiseProfile write SetNoiseProfile;
// these have defaults
property Sensitivity: Double read FSensitivity write SetSensitivity;
property Gain: Double read FNoiseGain write SetGain;
property FreqSmoothingHz: Double read FFreqSmoothingHz write SetFreqSmoothingHz;
property AttackDecayTime: Double read FAttackDecayTime write SetAttackDecayTime;
property LeaveNoise: Boolean read FbLeaveNoise write FbLeaveNoise;
// don't mess with these.
property SampleRate: Double read FSampleRate;// write FSampleRate;
property WindowSize: Integer read FWindowSize;// write FWindowSize;
property SpectrumSize: Integer read FSpectrumSize;// write FSpectrumSize;
property MinSignalTime: Single read FMinSignalTime;// write FMinSignalTime; // in secs
// This must be assigned or av's will occur
property WriteProc: TNoiseWriteProc read FWriteProc write FWriteProc;
end;
type
{ TNoiseRemovalChannel }
TNoiseRemovalChannel = class(TNoiseRemoval, IPAIODataIOInterface)
HasProfile: Boolean;
ProfileComplete: Boolean;
procedure WriteDataIO(ASender: IPAIODataIOInterface; AData: PSingle; ASamples: Integer);
end;
{ TNoiseRemovalMultiChannel }
TNoiseRemovalMultiChannel = class(IPAIODataIOInterface)
private
FChannels,
FSampleRate: Integer;
FHelper: TPAIOChannelHelper;
FNoise: array of TNoiseRemovalChannel;
FWriteProc: TNoiseWriteProc;
//IPAIODataIOInterface
procedure WriteDataIO(ASender: IPAIODataIOInterface; AData: PSingle; ASamples: Integer);
procedure DataWrite(ASender: TObject; AData: PSingle; ASampleCount: Integer);
public
constructor Create(AChannels: Integer; ASampleRate: Integer);
destructor Destroy; override;
procedure ReadNoiseProfile(AData: PSingle; ASamples: Integer);
procedure ProcessNoise(AData: PSingle; ASamples: Integer);
procedure Flush;
property WriteProc: TNoiseWriteProc read FWriteProc write FWriteProc;
end;
{ TuosNoiseRemoval }
type
TuosNoiseRemoval = class(TNoiseRemovalMultiChannel)
OutStream: TStream;
public
isprofiled : boolean ;
samprate : integer ;
procedure WriteData(ASender: TObject; AData: PSingle; ASampleCount: Integer) ;
function FilterNoise(ANoisyAudio: PSingle; InFrames: Integer; out Samples: Integer): PSingle;
end;
implementation
uses
math;
const
PI = 3.14159265358979323846;
MAX_HFFT = 10;
var
FFTArray: array[0..MAX_HFFT-1] of PFFT;
FFTLockCount: array[0..MAX_HFFT-1] of Integer;
{ TMultiChannelNoiseRemoval }
procedure TNoiseRemovalMultiChannel.WriteDataIO(ASender: IPAIODataIOInterface; AData: PSingle; ASamples: Integer);
begin
if Assigned(FWriteProc) then
FWriteProc(Self, AData, ASamples);
end;
procedure TNoiseRemovalMultiChannel.DataWrite(ASender: TObject; AData: PSingle; ASampleCount: Integer);
begin
(FHelper as IPAIODataIOInterface).WriteDataIO(ASender as IPAIODataIOInterface, AData, ASampleCount);
end;
constructor TNoiseRemovalMultiChannel.Create(AChannels: Integer;
ASampleRate: Integer);
var
i: Integer;
begin
FChannels:=AChannels;
FSampleRate:=ASampleRate;
FHelper := TPAIOChannelHelper.Create(Self);
SetLength(FNoise, AChannels);
for i := 0 to High(FNoise) do
begin
FNoise[i] := TNoiseRemovalChannel.Create;
FNoise[i].WriteProc:=@DataWrite;
FNoise[i].Init(ASampleRate);
FHelper.Outputs.Add(FNoise[i] as IPAIODataIOInterface);
end;
end;
destructor TNoiseRemovalMultiChannel.Destroy;
var
i: Integer;
begin
for i := 0 to High(FNoise) do
begin
FNoise[i].Free;
end;
SetLength(FNoise, 0);
FHelper.Free;
end;
procedure TNoiseRemovalMultiChannel.ReadNoiseProfile(AData: PSingle;
ASamples: Integer);
var
i: Integer;
begin
FHelper.Write(AData, ASamples);
for i := 0 to High(FNoise) do
begin
FNoise[i].ProfileComplete:=True;
FNoise[i].Process(nil, 0, True, False);
FNoise[i].HasProfile:=True;
FNoise[i].Init(FSampleRate);
end;
end;
procedure TNoiseRemovalMultiChannel.ProcessNoise(AData: PSingle;
ASamples: Integer);
begin
FHelper.Write(AData, ASamples);
end;
procedure TNoiseRemovalMultiChannel.Flush;
var
i: Integer;
begin
for i := 0 to High(FNoise) do
FNoise[i].Flush;
end;
procedure TNoiseRemovalChannel.WriteDataIO(ASender: IPAIODataIOInterface;
AData: PSingle; ASamples: Integer);
begin
Process(AData, ASamples, not HasProfile, not HasProfile);
end;
{ TuosNoiseRemoval }
procedure TuosNoiseRemoval.WriteData(ASender: TObject; AData: PSingle;
ASampleCount: Integer);
begin
OutStream.Write(AData^, ASampleCount*SizeOf(Single));
end;
function TuosNoiseRemoval.FilterNoise(ANoisyAudio: PSingle; InFrames: Integer; out Samples: Integer): PSingle;
var
MNoisyAudio: TMemoryStream;
begin
OutStream := TMemoryStream.Create;
MNoisyAudio := TMemoryStream.Create;
MNoisyAudio.Write(ANoisyAudio^, InFrames*SizeOf(Single));
MNoisyAudio.Position:=0;
if isprofiled = false then // take the first chunk as noisy sample
begin
ReadNoiseProfile(PSingle(MNoisyAudio.Memory), MNoisyAudio.Size div SizeOf(Single));
isprofiled := true;
end;
Result := nil;
ProcessNoise(PSingle(MNoisyAudio.Memory), MNoisyAudio.Size div SizeOf(Single));
Result:=GetMem(OutStream.Size);
Samples := OutStream.Size div SizeOf(Single);
OutStream.Position:=0;
OutStream.Read(Result^, OutStream.Size);
MNoisyAudio.free;
OutStream.Free;
// Result := ANoisyAudio;
end;
{ TNoiseRemoval }
function NewFloatArray(ALength: Integer): PSingle; inline;
begin
Result := Getmem(ALength*SizeOf(Single));
end;
procedure TNoiseRemoval.Initialize;
var
i: Integer;
begin
FFreqSmoothingBins := Trunc(FFreqSmoothingHz * FWindowSize / FSampleRate);
FAttackDecayBlocks := 1 + Trunc(FAttackDecayTime * FSampleRate / (FWindowSize / 2));
// Applies to amplitudes, divide by 20:
FNoiseAttenFactor := power(10, FNoiseGain/20);
// Applies to gain factors which apply to amplitudes, divide by 20:
//FOneBlockAttackDecay := power(10.0, (FNoiseGain / (20.0 * FAttackDecayBlocks)));
FOneBlockAttackDecay := power(10.0, (FNoiseGain / FAttackDecayBlocks) / 20 );
// Applies to power, divide by 10:
FSensitivityFactor := power(10.0, FSensitivity/10.0);
FMinSignalBlocks := Trunc(FMinSignalTime * FSampleRate / (FWindowSize / 2));
if( FMinSignalBlocks < 1 ) then
FMinSignalBlocks := 1;
FHistoryLen := (2 * FAttackDecayBlocks) - 1;
if (FHistoryLen < FMinSignalBlocks) then
FHistoryLen := FMinSignalBlocks;
SetLength(FSpectrums, FHistoryLen);
SetLength(FGains, FHistoryLen);
SetLength(FRealFFTs, FHistoryLen);
SetLength(FImagFFTs, FHistoryLen);
for i := 0 to FHistoryLen-1 do
begin
FSpectrums[i] := NewFloatArray(FSpectrumSize);
FGains[i] := NewFloatArray(FSpectrumSize);
FRealFFTs[i] := NewFloatArray(FSpectrumSize);
FImagFFTs[i] := NewFloatArray(FSpectrumSize);
end;
// Initialize the FFT
FFFT := TFFT.InitializeFFT(FWindowSize);
FFFTBuffer := NewFloatArray(FWindowSize);
FInWaveBuffer := NewFloatArray(FWindowSize);
FWindow := NewFloatArray(FWindowSize);
FOutOverlapBuffer := NewFloatArray(FWindowSize);
// Create a Hanning window function
for i := 0 to FWindowSize-1 do
FWindow[i] := 0.5 - 0.5 * cos((2.0*pi*i) / FWindowSize);
if FDoProfile then
begin
FillChar(FNoiseThreshold[0], SizeOf(FNoiseThreshold[0])*FSpectrumSize, 0);
//for i := 0 to FSpectrumSize-1 do
// FNoiseThreshold[i] := float(0);
end;
end;
procedure TNoiseRemoval.Reset;
var
i, j: Integer;
begin
for i := 0 to FHistoryLen-1 do
begin
for j := 0 to FSpectrumSize-1 do
begin
FSpectrums[i][j] := 0;
FGains[i][j] := FNoiseAttenFactor;
FRealFFTs[i][j] := 0.0;
FImagFFTs[i][j] := 0.0;
end;
end;
for j := 0 to FWindowSize-1 do
FOutOverlapBuffer[j] := 0.0;
FInputPos := 0;
FInSampleCount := 0;
FOutSampleCount := -(FWindowSize div 2) * (FHistoryLen - 1);
end;
function Min(A, B: Integer): Integer;
begin
if A < B then
Exit(A);
Result := B;
end;
procedure TNoiseRemoval.ProcessSamples(len: Integer; Buffer: PSingle);
var
i: Integer;
avail: Integer;
begin
//while((len and FOutSampleCount) < FInSampleCount) do
while len > 0 do
begin
avail := Min(len, FWindowSize - FInputPos);
for i := 0 to avail-1 do
FInWaveBuffer[FInputPos + i] := buffer[i];
buffer += avail;
len -= avail;
FInputPos += avail;
if (FInputPos = FWindowSize) then
begin
FillFirstHistoryWindow();
if (FDoProfile) then
GetProfile()
else
RemoveNoise();
RotateHistoryWindows();
// Rotate halfway for overlap-add
//for(i = 0; i < mWindowSize / 2; i++) {
for i := 0 to FWindowSize div 2 -1 do
FInWaveBuffer[i] := FInWaveBuffer[i + FWindowSize div 2];
FInputPos := FWindowSize div 2;
end;
end;
end;
procedure TNoiseRemoval.FillFirstHistoryWindow;
var
i: Integer;
begin
for i := 0 to FWindowSize-1 do
FFFTBuffer[i] := FInWaveBuffer[i];
FFFT^.RealFFTf(FFFTBuffer);
//for(i = 1; i < (mSpectrumSize-1); i++) {
for i := 1 to FSpectrumSize-2 do
begin
FRealFFTs[0][i] := FFFTBuffer[FFFT^.BitReversed[i] ];
FImagFFTs[0][i] := FFFTBuffer[FFFT^.BitReversed[i]+1];
FSpectrums[0][i] := FRealFFTs[0][i]*FRealFFTs[0][i] + FImagFFTs[0][i]*FImagFFTs[0][i];
FGains[0][i] := FNoiseAttenFactor;
end;
// DC and Fs/2 bins need to be handled specially
FSpectrums[0][0] := FFFTBuffer[0]*FFFTBuffer[0];
FSpectrums[0][FSpectrumSize-1] := FFFTBuffer[1]*FFFTBuffer[1];
FGains[0][0] := FNoiseAttenFactor;
FGains[0][FSpectrumSize-1] := FNoiseAttenFactor;
end;
function Max(A,B: Integer): Integer; inline;
begin
if A>B then
Exit(A);
Result := B;
end;
procedure TNoiseRemoval.ApplyFreqSmoothing(ASpec: PSingle);
var
tmp: PSingle;
i, j, j0, j1: Integer;
begin
tmp := NewFloatArray(FSpectrumSize);
for i := 0 to FSpectrumSize-1 do
begin
j0 := Max(0, i - FFreqSmoothingBins);
j1 := Min(FSpectrumSize-1, i + FFreqSmoothingBins);
tmp[i] := 0.0;
//for(j = j0; j <= j1; j++)
for j := j0 to j1-1 do
begin
tmp[i] += Aspec[j];
end;
tmp[i] := tmp[i] / (j1 - j0 + 1);
end;
//for(i = 0; i < mSpectrumSize; i++)
for i := 0 to FSpectrumSize-1 do
Aspec[i] := tmp[i];
Freemem(Tmp);
end;
procedure TNoiseRemoval.GetProfile;
var
start,
finish,
i, j: Integer;
min: Single;
begin
// The noise threshold for each frequency is the maximum
// level achieved at that frequency for a minimum of
// mMinSignalBlocks blocks in a row - the max of a min.
start := FHistoryLen - FMinSignalBlocks;
finish := FHistoryLen;
for j := 0 to FSpectrumSize-1 do
begin
min := FSpectrums[start][j];
for i := start+1 to finish-1 do
if (FSpectrums[i][j] < min) then
min := FSpectrums[i][j];
if (min > FNoiseThreshold[j]) then
FNoiseThreshold[j] := min;
end;
FOutSampleCount += FWindowSize div 2; // what is this for? Not used when we are getting the profile?
end;
procedure TNoiseRemoval.RemoveNoise;
var
center: Integer;
start,
finish,
i,j : Integer;
min: Single;
out_: Integer;
begin
center := FHistoryLen div 2;
start := center - FMinSignalBlocks div 2;
finish := start + FMinSignalBlocks;
// Raise the gain for elements in the center of the sliding history
for j := 0 to FSpectrumSize-1 do
begin
min := FSpectrums[start][j];
//for (i = start+1; i < finish; i++) {
for i := start+1 to finish-1 do
begin
if (FSpectrums[i][j] < min) then
min := FSpectrums[i][j];
end;
if (min > FSensitivityFactor * FNoiseThreshold[j]) and (FGains[center][j] < 1.0) then
begin
if (FbLeaveNoise) then
FGains[center][j] := 0.0
else
FGains[center][j] := 1.0;
end
else
begin
if (FbLeaveNoise) then
FGains[center][j] := 1.0;
end;
end;
// Decay the gain in both directions;
// note that mOneBlockAttackDecay is less than 1.0
// of linear attenuation per block
for j := 0 to FSpectrumSize-1 do
begin
for i := center+1 to FHistoryLen-1 do
begin
if (FGains[i][j] < FGains[i - 1][j] * FOneBlockAttackDecay) then
FGains[i][j] := FGains[i - 1][j] * FOneBlockAttackDecay;
if (FGains[i][j] < FNoiseAttenFactor) then
FGains[i][j] := FNoiseAttenFactor;
end;
for i := center-1 downto 0 do
begin
if (FGains[i][j] < FGains[i + 1][j] * FOneBlockAttackDecay) then
FGains[i][j] := FGains[i + 1][j] * FOneBlockAttackDecay;
if (FGains[i][j] < FNoiseAttenFactor) then
FGains[i][j] := FNoiseAttenFactor;
end;
end;
// Apply frequency smoothing to output gain
out_ := FHistoryLen - 1; // end of the queue
ApplyFreqSmoothing(FGains[out_]);
// Apply gain to FFT
//for (j = 0; j < (mSpectrumSize-1); j++) {
for j := 0 to FSpectrumSize-2 do
begin
FFFTBuffer[j*2 ] := FRealFFTs[out_][j] * FGains[out_][j];
FFFTBuffer[j*2+1] := FImagFFTs[out_][j] * FGains[out_][j];
end;
// The Fs/2 component is stored as the imaginary part of the DC component
FFFTBuffer[1] := FRealFFTs[out_][FSpectrumSize-1] * FGains[out_][FSpectrumSize-1];
// Invert the FFT into the output buffer
FFFT^.InverseRealFFTf(FFFTBuffer);
// Overlap-add
for j := 0 to FSpectrumSize-2 do
begin
FOutOverlapBuffer[j*2 ] += FFFTBuffer[FFFT^.BitReversed[j] ] * FWindow[j*2 ];
FOutOverlapBuffer[j*2+1] += FFFTBuffer[FFFT^.BitReversed[j]+1] * FWindow[j*2+1];
end;
// Output the first half of the overlap buffer, they're done -
// and then shift the next half over.
if (FOutSampleCount >= 0) then // ...but not if it's the first half-window
begin
//FOutputTrack->Append((samplePtr)mOutOverlapBuffer, floatSample, mWindowSize / 2);
FWriteProc(Self, FOutOverlapBuffer, FWindowSize div 2);
end;
FOutSampleCount += FWindowSize div 2;
//for(j = 0; j < mWindowSize / 2; j++)
for j := 0 to FWindowSize div 2 -1 do
begin
FOutOverlapBuffer[j] := FOutOverlapBuffer[j + (FWindowSize div 2)];
FOutOverlapBuffer[j + (FWindowSize div 2)] := 0.0;
end
end;
procedure TNoiseRemoval.RotateHistoryWindows;
var
last: Integer;
i: Integer;
lastSpectrum: PSingle;
lastGain: PSingle;
lastRealFFT: PSingle;
lastImagFFT: PSingle;
begin
last := FHistoryLen - 1;
// Remember the last window so we can reuse it
lastSpectrum := FSpectrums[last];
lastGain := FGains[last];
lastRealFFT := FRealFFTs[last];
lastImagFFT := FImagFFTs[last];
// Rotate each window forward
//for(i = last; i >= 1; i--) {
for i := last downto 1 do
begin
FSpectrums[i] := FSpectrums[i-1];
FGains[i] := FGains[i-1];
FRealFFTs[i] := FRealFFTs[i-1];
FImagFFTs[i] := FImagFFTs[i-1];
end;
// Reuse the last buffers as the new first window
FSpectrums[0] := lastSpectrum;
FGains[0] := lastGain;
FRealFFTs[0] := lastRealFFT;
FImagFFTs[0] := lastImagFFT;
end;
procedure TNoiseRemoval.FinishTrack;
var
empty: PSingle;
i: Integer;
begin
// Keep flushing empty input buffers through the history
// windows until we've output exactly as many samples as
// were input.
// Well, not exactly, but not more than mWindowSize/2 extra samples at the end.
// We'll delete them later in ProcessOne.
empty := NewFloatArray(FWindowSize div 2);
//for(i = 0; i < mWindowSize / 2; i++)
for i := 0 to FWindowSize div 2 -1 do
empty[i] := 0.0;
while (FOutSampleCount < FInSampleCount) do
ProcessSamples(FWindowSize div 2, empty);
Freemem(empty);
end;
procedure TNoiseRemoval.Cleanup;
var
i: Integer;
begin
FFFT^.EndFFT;
if (FDoProfile) then
ApplyFreqSmoothing(@FNoiseThreshold[0]);
for i := 0 to FHistoryLen-1 do
begin
FreeMem(FSpectrums[i]);
FreeMem(FGains[i]);
FreeMem(FRealFFTs[i]);
FreeMem(FImagFFTs[i]);
end;
SetLength(FSpectrums,0);
SetLength(FGains,0);
SetLength(FRealFFTs,0);
SetLength(FImagFFTs,0);
FreeMem(FFFTBuffer);
FreeMem(FInWaveBuffer);
FreeMem(FWindow);
FreeMem(FOutOverlapBuffer);
FInited := False;
end;
function TNoiseRemoval.GetNoiseProfile: TSingleArray;
begin
SetLength(Result, FSpectrumSize);
Move(FNoiseThreshold[0], Result[0], FSpectrumSize);
end;
procedure TNoiseRemoval.SetAttackDecayTime(AValue: Double);
begin
if FAttackDecayTime=AValue then Exit;
if AValue < 0.0 then AValue := 0;
if AValue > 1.0 then AValue := 1.0;
FAttackDecayTime:=AValue;
end;
procedure TNoiseRemoval.SetFreqSmoothingHz(AValue: Double);
begin
if FFreqSmoothingHz=AValue then Exit;
if AValue<0 then AValue:=0;
if AValue>1000 then AValue := 1000;
FFreqSmoothingHz:=AValue;
end;
procedure TNoiseRemoval.SetGain(AValue: Double);
begin
if FNoiseGain=AValue then Exit;
if AValue > 0 then AValue:=0;
if AValue < -48 then AValue := -48;
FNoiseGain:=AValue;
end;
procedure TNoiseRemoval.SetNoiseProfile(AValue: TSingleArray);
begin
SetLength(FNoiseThreshold, FSpectrumSize);
Move(AValue[0], FNoiseThreshold[0], FSpectrumSize);
FHasProfile:=True;
FDoProfile:=False;
Cleanup; // set after FDoProfile so the profile is not processed.
end;
procedure TNoiseRemoval.SetSensitivity(AValue: Double);
begin
if FSensitivity=AValue then Exit;
if AValue < -20.0 then AValue:=-20.0;
if AValue > 20.0 then AValue := 20.0;
FSensitivity:=AValue;
end;
constructor TNoiseRemoval.Create;
begin
FWindowSize:=2048;
FSpectrumSize:= 1 + FWindowSize div 2;
// loaded prefs
FSensitivity := 0.0;
FNoiseGain := -24.0;
FFreqSmoothingHz := 150.0;
FAttackDecayTime:= 0.15;
FbLeaveNoise:=False;
FMinSignalTime := 0.05;
FHasProfile := False;
FDoProfile := True;
SetLength(FNoiseThreshold, FSpectrumSize);
end;
destructor TNoiseRemoval.Destroy;
begin
SetLength(FNoiseThreshold, 0);
inherited Destroy;
end;
function TNoiseRemoval.Init(ASampleRate: Integer): Boolean;
begin
FSampleRate:=ASampleRate;
Initialize;
FInited:=True;
Result := True;
Reset;
end;
function TNoiseRemoval.Process(AData: PSingle; ASampleCount: Integer;
AGetNoiseProfile: Boolean; AMoreComing: Boolean = False): Boolean;
begin
if not FInited then
Raise Exception.Create('TNoiseRemoval is not Inited');
if not AGetNoiseProfile and not FHasProfile then
raise Exception.Create('Tried to remove noise without profile.');
FDoProfile:=AGetNoiseProfile;
if FDoProfile and (FTotalRead = 0) then
begin
Initialize;
reset;
end;
Inc(FTotalRead, ASampleCount);
ProcessSamples(ASampleCount, AData);
Result := True;
if AMoreComing then
Exit;
if {AGetNoiseProfile or }FDoProfile then
begin
FHasProfile:=True;
Cleanup; // triggers the data in FNoiseThreshold to be processed
// must be set after Cleanup() is called
//FDoProfile:=False;
//Initialize;
FTotalRead:=0;
end;
FHasProfile:=True;
FDoProfile := False;
end;
procedure TNoiseRemoval.Flush;
begin
if not FInited then
Exit;
FinishTrack;
Cleanup; // Sets FInited to False
end;
{ TFFT }
function TFFT.InitializeFFT(FFTLen: Integer): PFFT;
var
i: Integer;
temp: Integer;
mask: Integer;
begin
Result := New(PFFT);
if Result = nil then
raise EOutOfMemory.Create('Error allocating memory for FFT');
with Result^ do begin
{*
* FFT size is only half the number of data points
* The full FFT output can be reconstructed from this FFT's output.
* (This optimization can be made since the data is real.)
*}
Points := FFTLen div 2;
SetLength(FPCSinTable, 2*Points);
SinTable:=@FPCSinTable[0];
BitReversed := Getmemory(Points*SizeOf(Integer));
if BitReversed = nil then
raise EOutOfMemory.Create('Error allocating memory for BitReversed.');
for i := 0 to Points-1 do
begin
temp:=0;
mask := Points div 2;
while mask > 0 do
begin
//for(mask=h->Points/2;mask>0;mask >>= 1)
// temp=(temp >> 1) + (i&mask ? h->Points : 0);
temp := (temp shr 1);
if (i and mask) <> 0 then
temp := temp + Points;
//else temp := temp + 0; // why would you do that?
mask := mask shr 1;
end;
BitReversed[i]:=temp;
end;
for i := 0 to Points-1 do
begin
SinTable[BitReversed[i] ]:= -sin(2*PI*i/(2*Points));
SinTable[BitReversed[i]+1]:= -cos(2*PI*i/(2*Points));
end;
{$ifdef EXPERIMENTAL_EQ_SSE_THREADED}
// new SSE FFT routines work on live data
for(i=0;i<32;i++)
if((1<<i)&fftlen)
h->pow2Bits=i;
{$endif}
end; // with Result^
end;
procedure TFFT.EndFFT;
begin
if Points>0 then
begin
Freemem(BitReversed);
SetLength(FPCSinTable, 0);
SinTable:=nil;
end;
Dispose(PFFT(@Self));
end;
function TFFT.GetFFT(FFTLen: Integer): PFFT;
var
h: Integer = 0;
n: Integer;
begin
n := fftlen div 2;
while (h<MAX_HFFT) and (FFTArray[h] <> nil) and (n <> FFTArray[h]^.Points) do
begin
if (h<MAX_HFFT) then
begin
if(FFTArray[h] = nil) then
begin
FFTArray[h] := InitializeFFT(fftlen);
FFTLockCount[h] := 0;
end;
Inc(FFTLockCount[h]);
Exit(FFTArray[h]);
end
else begin
// All buffers used, so fall back to allocating a new set of tables
Exit(InitializeFFT(fftlen));
end;
Inc(h);
end;
end;
procedure TFFT.ReleaseFFT;
var
h: Integer = 0;
begin
while (h<MAX_HFFT) and (FFTArray[h] <> @Self) do
begin
if(h<MAX_HFFT) then
begin
Dec(FFTLockCount[h]);
end
else
begin
EndFFT;
end;
Inc(h);
end;
end;
procedure TFFT.InverseRealFFTf(buffer: PSingle);
var
A, B: PSingle;
sptr: PSingle;
endptr1, endptr2: PSingle;
br1: PInteger;
HRplus,HRminus,HIplus,HIminus: Single;
v1,v2,sin,cos: Single;
ButterfliesPerGroup: Integer;
begin
ButterfliesPerGroup:=Points div 2;
//* Massage input to get the input for a real output sequence. */
A:=@buffer[2];
B:=@buffer[Points*2-2];
br1:=@BitReversed[1];
while(A<B) do
begin
sin:=SinTable[br1^];
cos:=SinTable[br1[1]];
//HRplus = (HRminus = *A - *B ) + (*B * 2);
HRminus:=A^-B^;
HRplus:=HRminus+ (B^ *2);
//HIplus = (HIminus = *(A+1) - *(B+1)) + (*(B+1) * 2);
HIminus:=A[1]-B[1];
HIplus:=HIminus+(B[1] *2);
v1 := (sin*HRminus + cos*HIplus);
v2 := (cos*HRminus - sin*HIplus);
A^ := (HRplus + v1) * single(0.5);
B^ := A^ - v1;
A[1] := (HIminus - v2) * single(0.5);
B[1] := A[1] - HIminus;
A+=2;
B-=2;
Inc(br1);
end;
//* Handle center bin (just need conjugate) */
A[1] :=-A[1];
{* Handle DC bin separately - this ignores any Fs/2 component
buffer[1]=buffer[0]=buffer[0]/2;*}
//* Handle DC and Fs/2 bins specially */
//* The DC bin is passed in as the real part of the DC complex value */
//* The Fs/2 bin is passed in as the imaginary part of the DC complex value */
//* (v1+v2) = buffer[0] == the DC component */
//* (v1-v2) = buffer[1] == the Fs/2 component */
v1:=0.5*(buffer[0]+buffer[1]);
v2:=0.5*(buffer[0]-buffer[1]);
buffer[0]:=v1;
buffer[1]:=v2;
{*
* Butterfly:
* Ain-----Aout
* \ /
* / \
* Bin-----Bout
*}
endptr1:=@buffer[Points*2];
while(ButterfliesPerGroup>0) do
begin
A:=buffer;
B:=@buffer[ButterfliesPerGroup*2];
sptr:=@SinTable[0];
while(A<endptr1) do
begin
sin:=sptr^; Inc(sptr); // *(sptr++);
cos:=sptr^; Inc(sptr); // *(sptr++);
endptr2:=B;
while(A<endptr2) do
begin
v1:=B^*cos - B[1]*sin;
v2:=B^*sin + B[1]*cos;
B^ := (A^+v1)*Single(0.5);
A^ := B^ - v1; Inc(A); Inc(B); //*(A++)=*(B++)-v1;
B^ := (A^ + v2)* Single(0.5); //*B=(*A+v2)*(fft_type)0.5;
A^ := B^ - v2; Inc(A); Inc(B); //*(A++)=*(B++)-v2;
end;
A:=B;
B := @B[ButterfliesPerGroup*2];
end;
ButterfliesPerGroup := ButterfliesPerGroup shr 1;
end;
end;
procedure TFFT.CleanupFFT;
var
h: Integer;
begin
for h :=0 to MAX_HFFT-1do begin
if((FFTLockCount[h] <= 0) and (FFTArray[h] <> nil)) then
begin
FFTArray[h]^.EndFFT;
FFTArray[h] := nil;
end;
end;
end;
procedure TFFT.RealFFTf(buffer: PSingle);
var
A, B: PSingle;
sptr: PSingle;
endptr1, endptr2: PSingle;
br1, br2: PInteger;
HRplus,HRminus,HIplus,HIminus: Single;
v1,v2,sin_,cos_: Single;
ButterfliesPerGroup: Integer;
begin
ButterfliesPerGroup:=Points div 2;
{*
* Butterfly:
* Ain-----Aout
* \ /
* / \
* Bin-----Bout
*}
endptr1:=buffer+Points*2;
while(ButterfliesPerGroup>0) do
begin
A:=buffer;
B:=buffer+ButterfliesPerGroup*2;
sptr:=@SinTable[0];
while(A<endptr1) do
begin
sin_:=sptr^;
cos_ := sptr[1];
endptr2:=B;
while(A<endptr2) do
begin
v1 := B^ * cos_ + B[1] * sin_; //v1=*B*cos + *(B+1)*sin;
v2 := B^ * sin_ - B[1] * cos_; //v2=*B*sin - *(B+1)*cos;
B^ := A^+v1; //*B=(*A+v1);
A^ := B^-2*v1; Inc(A); Inc(B); //*(A++)=*(B++)-2*v1;
B^ := A^-v2; //*B=(*A-v2);
A^ := B^+2*v2; Inc(A); Inc(B); //*(A++)=*(B++)+2*v2;
end;
A:=B;
B:=B+ButterfliesPerGroup*2;
sptr:=sptr+2;
end;
ButterfliesPerGroup := ButterfliesPerGroup shr 1;
end;
//* Massage output to get the output for a real input sequence. */
br1:=@BitReversed[1]; // is this wrong? Should be @BitReversed[0] ; ?
br2:=@BitReversed[Points-1];
while(br1<br2) do
begin
sin_:=SinTable[br1[0]];
cos_:=SinTable[br1[1]];
A:=@buffer[br1^];
B:=@buffer[br2^];
//HRplus = (HRminus = *A - *B ) + (*B * 2);
HRminus := A^ - B^;
HRplus := HRminus + (B^ * 2);
//HIplus = (HIminus = *(A+1) - *(B+1)) + (*(B+1) * 2);
HIminus := A[1] - B[1];
HIplus := HIminus + (B[1] * 2);
v1 := (sin_*HRminus - cos_*HIplus);
v2 := (cos_*HRminus + sin_*HIplus);
A^ := (HRplus + v1) * single(0.5);
B^ := A^ - v1;
A[1] := (HIminus + v2) * single(0.5);
B[1] := A[1] - HIminus;
Inc(br1);
Dec(br2);
end;
//* Handle the center bin (just need a conjugate) */
A:=buffer+br1[1];
A^ := -A^;
{* Handle DC bin separately - and ignore the Fs/2 bin
buffer[0]+=buffer[1];
buffer[1]=(fft_type)0;*}
///* Handle DC and Fs/2 bins separately */
///* Put the Fs/2 value into the imaginary part of the DC bin */
v1:=buffer[0]-buffer[1];
buffer[0]+=buffer[1];
buffer[1]:=v1;
end;
procedure TFFT.ReorderToTime(Buffer: PSingle; TimeOut: PSingle);
var
i: Integer;
begin
// Copy the data into the real outputs
//for(int i=0;i<hFFT->Points;i++) {
for i := 0 to Points-1 do
begin
TimeOut[i*2 ]:=buffer[BitReversed[i] ];
TimeOut[i*2+1]:=buffer[BitReversed[i]+1];
end;
end;
procedure TFFT.ReorderToFreq(Buffer: PSingle; RealOut: PSingle; ImagOut: PSingle
);
var
i: Integer;
begin
// Copy the data into the real and imaginary outputs
//for(int i=1;i<hFFT->Points;i++)
for i := 1 to Points-1 do
begin
RealOut[i]:=buffer[BitReversed[i] ];
ImagOut[i]:=buffer[BitReversed[i]+1];
end;
RealOut[0] := buffer[0]; // DC component
ImagOut[0] := 0;
RealOut[Points] := buffer[1]; // Fs/2 component
ImagOut[Points] := 0;
end;
initialization
FillChar(FFTArray, SizeOf(FFTArray), 0);
FillChar(FFTLockCount, SizeOf(FFTLockCount), 0);
end.
Reference in New Issue View Git Blame Copy Permalink