source: trunk/Packages/uos/uos_dsp_noiseremoval.pas

Last change on this file was 664, checked in by chronos, 3 days ago
  • Added: Ability to play music in background in start screen and in-game. Used uos as audio library.
  • Property svn:executable set to *
File size: 31.0 KB
Line 
1{This unit is part of United Openlibraries of Sound (uos)}
2{
3 audacity which this file was converted from is GPLv2
4 http://www.audacityteam.org/about/license/
5 Converted By Andrew Haines
6 License : modified LGPL.
7 Fred van Stappen fiens@hotmail.com
8 }
9
10unit uos_dsp_noiseremoval;
11
12{$mode objfpc}{$H+}
13
14interface
15
16uses
17 Classes, SysUtils, uos_dsp_utils ;
18
19type
20 PFFT = ^TFFT;
21
22 { TFFT }
23
24 TFFT = object
25 BitReversed: PInteger;
26 SinTable: PSingle;
27 Points: Integer;
28 FPCSinTable: array of Single;
29 function InitializeFFT(FFTLen: Integer): PFFT; static;
30 procedure EndFFT;
31 function GetFFT(FFTLen: Integer): PFFT; static;
32 procedure ReleaseFFT;
33 procedure InverseRealFFTf(buffer: PSingle);
34 procedure CleanupFFT; static; // ???
35 procedure RealFFTf(buffer: PSingle);
36 procedure ReorderToTime(Buffer: PSingle; TimeOut: PSingle);
37 procedure ReorderToFreq(Buffer: PSingle; RealOut: PSingle; ImagOut: PSingle);
38 end;
39
40type
41 PPSingle=^PSingle;
42 TSingleArray = array of Single;
43
44 TNoiseRemoval = class;
45 TNoiseWriteProc = procedure(ASender: TObject; AData: PSingle; ASampleCount: Integer) of Object;
46
47 { TNoiseRemoval }
48
49 TNoiseRemoval = class
50 private
51 FDoProfile: Boolean;
52 FHasProfile: Boolean;
53
54 // Parameters chosen before the first phase
55 FSampleRate: Double;
56 FWindowSize: Integer;
57 FSpectrumSize: Integer;
58 FMinSignalTime: Single; // in secs
59
60 // The frequency-indexed noise threshold derived during the first
61 // phase of analysis
62 FNoiseThreshold: array of Single; // length in FSpectrumSize
63
64 // Parameters that affect the noise removal, regardless of how the
65 // noise profile was extracted
66 FSensitivity: Double;
67 FFreqSmoothingHz: Double;
68 FNoiseGain: Double; // in dB, should be negative
69 FAttackDecayTime: Double;
70 FbLeaveNoise: Boolean;
71
72 // change this later
73 procedure Initialize;
74 procedure Reset; // StartNewTrack
75 procedure ProcessSamples(len: Integer; Buffer:PSingle);
76 procedure FillFirstHistoryWindow;
77 procedure ApplyFreqSmoothing(ASpec: PSingle);
78 procedure GetProfile;
79 procedure RemoveNoise;
80 procedure RotateHistoryWindows;
81 procedure FinishTrack;
82 procedure Cleanup;
83 private
84 //FOutputTrack: PSingle; // WaveTrack;
85 FInSampleCount: Integer;
86 FOutSampleCount: Integer;
87 FInputPos: Integer;
88 FFFT: PFFT;
89 FFFTBuffer: PSingle; // FWindowSize
90 FWindow: PSingle; // FWindowSize
91 FFreqSmoothingBins: Integer;
92 FAttackDecayBlocks: Integer;
93 FOneBlockAttackDecay: Single;
94 FNoiseAttenFactor: Single;
95 FSensitivityFactor: Single;
96 FMinSignalBlocks: Integer;
97 FHistoryLen: Integer;
98 FInWaveBuffer: PSingle; // FWindowSize
99 FOutOverlapBuffer: PSingle; // FWindowSize
100 FSpectrums: array of PSingle; // FHistoryLen x FSpectrumSize
101 FGains: array of PSingle; // FHistoryLen x FSpectrumSize
102 FRealFFTs: array of PSingle; // FHistoryLen x FWindowSize
103 FImagFFTs: array of PSingle; // FHistoryLen x FWindowSize
104 FWriteProc: TNoiseWriteProc;
105 FInited: Boolean;
106 FTotalRead: QWord;
107 function GetNoiseProfile: TSingleArray;
108 procedure SetAttackDecayTime(AValue: Double);
109 procedure SetFreqSmoothingHz(AValue: Double);
110 procedure SetGain(AValue: Double);
111 procedure SetNoiseProfile(AValue: TSingleArray);
112 procedure SetSensitivity(AValue: Double);
113 public
114 constructor Create;
115 destructor Destroy; override;
116 function Init(ASampleRate: Integer): Boolean;
117 function Process(AData: PSingle; ASampleCount: Integer; AGetNoiseProfile: Boolean; AMoreComing: Boolean = False): Boolean;
118 procedure Flush; // finish writing out data in buffers.
119
120 property NoiseProfile: TSingleArray read GetNoiseProfile write SetNoiseProfile;
121
122 // these have defaults
123 property Sensitivity: Double read FSensitivity write SetSensitivity;
124 property Gain: Double read FNoiseGain write SetGain;
125 property FreqSmoothingHz: Double read FFreqSmoothingHz write SetFreqSmoothingHz;
126 property AttackDecayTime: Double read FAttackDecayTime write SetAttackDecayTime;
127 property LeaveNoise: Boolean read FbLeaveNoise write FbLeaveNoise;
128
129 // don't mess with these.
130 property SampleRate: Double read FSampleRate;// write FSampleRate;
131 property WindowSize: Integer read FWindowSize;// write FWindowSize;
132 property SpectrumSize: Integer read FSpectrumSize;// write FSpectrumSize;
133 property MinSignalTime: Single read FMinSignalTime;// write FMinSignalTime; // in secs
134
135 // This must be assigned or av's will occur
136 property WriteProc: TNoiseWriteProc read FWriteProc write FWriteProc;
137 end;
138
139type
140
141 { TNoiseRemovalChannel }
142
143 TNoiseRemovalChannel = class(TNoiseRemoval, IPAIODataIOInterface)
144 HasProfile: Boolean;
145 ProfileComplete: Boolean;
146 procedure WriteDataIO(ASender: IPAIODataIOInterface; AData: PSingle; ASamples: Integer);
147 end;
148
149 { TNoiseRemovalMultiChannel }
150
151 TNoiseRemovalMultiChannel = class(IPAIODataIOInterface)
152 private
153 FChannels,
154 FSampleRate: Integer;
155 FHelper: TPAIOChannelHelper;
156 FNoise: array of TNoiseRemovalChannel;
157 FWriteProc: TNoiseWriteProc;
158 //IPAIODataIOInterface
159 procedure WriteDataIO(ASender: IPAIODataIOInterface; AData: PSingle; ASamples: Integer);
160 procedure DataWrite(ASender: TObject; AData: PSingle; ASampleCount: Integer);
161 public
162 constructor Create(AChannels: Integer; ASampleRate: Integer);
163 destructor Destroy; override;
164 procedure ReadNoiseProfile(AData: PSingle; ASamples: Integer);
165 procedure ProcessNoise(AData: PSingle; ASamples: Integer);
166 procedure Flush;
167 property WriteProc: TNoiseWriteProc read FWriteProc write FWriteProc;
168 end;
169
170{ TuosNoiseRemoval }
171 type
172 TuosNoiseRemoval = class(TNoiseRemovalMultiChannel)
173 OutStream: TStream;
174 public
175 isprofiled : boolean ;
176 samprate : integer ;
177 procedure WriteData(ASender: TObject; AData: PSingle; ASampleCount: Integer) ;
178 function FilterNoise(ANoisyAudio: PSingle; InFrames: Integer; out Samples: Integer): PSingle;
179 end;
180
181 implementation
182uses
183 math;
184
185const
186 PI = 3.14159265358979323846;
187 MAX_HFFT = 10;
188var
189 FFTArray: array[0..MAX_HFFT-1] of PFFT;
190 FFTLockCount: array[0..MAX_HFFT-1] of Integer;
191
192{ TMultiChannelNoiseRemoval }
193
194procedure TNoiseRemovalMultiChannel.WriteDataIO(ASender: IPAIODataIOInterface; AData: PSingle; ASamples: Integer);
195begin
196 if Assigned(FWriteProc) then
197 FWriteProc(Self, AData, ASamples);
198end;
199
200procedure TNoiseRemovalMultiChannel.DataWrite(ASender: TObject; AData: PSingle; ASampleCount: Integer);
201begin
202 (FHelper as IPAIODataIOInterface).WriteDataIO(ASender as IPAIODataIOInterface, AData, ASampleCount);
203end;
204
205constructor TNoiseRemovalMultiChannel.Create(AChannels: Integer;
206 ASampleRate: Integer);
207var
208 i: Integer;
209begin
210 FChannels:=AChannels;
211 FSampleRate:=ASampleRate;
212 FHelper := TPAIOChannelHelper.Create(Self);
213 SetLength(FNoise, AChannels);
214 for i := 0 to High(FNoise) do
215 begin
216 FNoise[i] := TNoiseRemovalChannel.Create;
217 FNoise[i].WriteProc:=@DataWrite;
218 FNoise[i].Init(ASampleRate);
219 FHelper.Outputs.Add(FNoise[i] as IPAIODataIOInterface);
220 end;
221end;
222
223destructor TNoiseRemovalMultiChannel.Destroy;
224var
225 i: Integer;
226begin
227 for i := 0 to High(FNoise) do
228 begin
229 FNoise[i].Free;
230 end;
231 SetLength(FNoise, 0);
232 FHelper.Free;
233end;
234
235procedure TNoiseRemovalMultiChannel.ReadNoiseProfile(AData: PSingle;
236 ASamples: Integer);
237var
238 i: Integer;
239begin
240 FHelper.Write(AData, ASamples);
241 for i := 0 to High(FNoise) do
242 begin
243 FNoise[i].ProfileComplete:=True;
244 FNoise[i].Process(nil, 0, True, False);
245 FNoise[i].HasProfile:=True;
246 FNoise[i].Init(FSampleRate);
247 end;
248end;
249
250procedure TNoiseRemovalMultiChannel.ProcessNoise(AData: PSingle;
251 ASamples: Integer);
252begin
253 FHelper.Write(AData, ASamples);
254end;
255
256procedure TNoiseRemovalMultiChannel.Flush;
257var
258 i: Integer;
259begin
260 for i := 0 to High(FNoise) do
261 FNoise[i].Flush;
262end;
263
264procedure TNoiseRemovalChannel.WriteDataIO(ASender: IPAIODataIOInterface;
265 AData: PSingle; ASamples: Integer);
266begin
267 Process(AData, ASamples, not HasProfile, not HasProfile);
268end;
269
270{ TuosNoiseRemoval }
271
272procedure TuosNoiseRemoval.WriteData(ASender: TObject; AData: PSingle;
273 ASampleCount: Integer);
274begin
275 OutStream.Write(AData^, ASampleCount*SizeOf(Single));
276end;
277
278function TuosNoiseRemoval.FilterNoise(ANoisyAudio: PSingle; InFrames: Integer; out Samples: Integer): PSingle;
279var
280 MNoisyAudio: TMemoryStream;
281 begin
282 OutStream := TMemoryStream.Create;
283
284 MNoisyAudio := TMemoryStream.Create;
285 MNoisyAudio.Write(ANoisyAudio^, InFrames*SizeOf(Single));
286 MNoisyAudio.Position:=0;
287
288 if isprofiled = false then // take the first chunk as noisy sample
289 begin
290 ReadNoiseProfile(PSingle(MNoisyAudio.Memory), MNoisyAudio.Size div SizeOf(Single));
291 isprofiled := true;
292 end;
293
294 Result := nil;
295
296 ProcessNoise(PSingle(MNoisyAudio.Memory), MNoisyAudio.Size div SizeOf(Single));
297
298 Result:=GetMem(OutStream.Size);
299 Samples := OutStream.Size div SizeOf(Single);
300 OutStream.Position:=0;
301 OutStream.Read(Result^, OutStream.Size);
302
303 MNoisyAudio.free;
304 OutStream.Free;
305
306 // Result := ANoisyAudio;
307 end;
308
309{ TNoiseRemoval }
310
311function NewFloatArray(ALength: Integer): PSingle; inline;
312begin
313 Result := Getmem(ALength*SizeOf(Single));
314end;
315
316procedure TNoiseRemoval.Initialize;
317var
318 i: Integer;
319begin
320 FFreqSmoothingBins := Trunc(FFreqSmoothingHz * FWindowSize / FSampleRate);
321 FAttackDecayBlocks := 1 + Trunc(FAttackDecayTime * FSampleRate / (FWindowSize / 2));
322 // Applies to amplitudes, divide by 20:
323 FNoiseAttenFactor := power(10, FNoiseGain/20);
324
325 // Applies to gain factors which apply to amplitudes, divide by 20:
326 //FOneBlockAttackDecay := power(10.0, (FNoiseGain / (20.0 * FAttackDecayBlocks)));
327 FOneBlockAttackDecay := power(10.0, (FNoiseGain / FAttackDecayBlocks) / 20 );
328 // Applies to power, divide by 10:
329 FSensitivityFactor := power(10.0, FSensitivity/10.0);
330 FMinSignalBlocks := Trunc(FMinSignalTime * FSampleRate / (FWindowSize / 2));
331 if( FMinSignalBlocks < 1 ) then
332 FMinSignalBlocks := 1;
333 FHistoryLen := (2 * FAttackDecayBlocks) - 1;
334
335 if (FHistoryLen < FMinSignalBlocks) then
336 FHistoryLen := FMinSignalBlocks;
337
338 SetLength(FSpectrums, FHistoryLen);
339 SetLength(FGains, FHistoryLen);
340 SetLength(FRealFFTs, FHistoryLen);
341 SetLength(FImagFFTs, FHistoryLen);
342 for i := 0 to FHistoryLen-1 do
343 begin
344 FSpectrums[i] := NewFloatArray(FSpectrumSize);
345 FGains[i] := NewFloatArray(FSpectrumSize);
346 FRealFFTs[i] := NewFloatArray(FSpectrumSize);
347 FImagFFTs[i] := NewFloatArray(FSpectrumSize);
348 end;
349
350 // Initialize the FFT
351 FFFT := TFFT.InitializeFFT(FWindowSize);
352
353 FFFTBuffer := NewFloatArray(FWindowSize);
354 FInWaveBuffer := NewFloatArray(FWindowSize);
355 FWindow := NewFloatArray(FWindowSize);
356 FOutOverlapBuffer := NewFloatArray(FWindowSize);
357
358 // Create a Hanning window function
359 for i := 0 to FWindowSize-1 do
360 FWindow[i] := 0.5 - 0.5 * cos((2.0*pi*i) / FWindowSize);
361
362 if FDoProfile then
363 begin
364 FillChar(FNoiseThreshold[0], SizeOf(FNoiseThreshold[0])*FSpectrumSize, 0);
365 //for i := 0 to FSpectrumSize-1 do
366 // FNoiseThreshold[i] := float(0);
367 end;
368
369end;
370
371procedure TNoiseRemoval.Reset;
372var
373 i, j: Integer;
374begin
375 for i := 0 to FHistoryLen-1 do
376 begin
377 for j := 0 to FSpectrumSize-1 do
378 begin
379 FSpectrums[i][j] := 0;
380 FGains[i][j] := FNoiseAttenFactor;
381 FRealFFTs[i][j] := 0.0;
382 FImagFFTs[i][j] := 0.0;
383 end;
384 end;
385
386 for j := 0 to FWindowSize-1 do
387 FOutOverlapBuffer[j] := 0.0;
388
389 FInputPos := 0;
390 FInSampleCount := 0;
391 FOutSampleCount := -(FWindowSize div 2) * (FHistoryLen - 1);
392end;
393
394function Min(A, B: Integer): Integer;
395begin
396 if A < B then
397 Exit(A);
398 Result := B;
399end;
400
401procedure TNoiseRemoval.ProcessSamples(len: Integer; Buffer: PSingle);
402var
403 i: Integer;
404 avail: Integer;
405begin
406 //while((len and FOutSampleCount) < FInSampleCount) do
407 while len > 0 do
408 begin
409 avail := Min(len, FWindowSize - FInputPos);
410 for i := 0 to avail-1 do
411 FInWaveBuffer[FInputPos + i] := buffer[i];
412 buffer += avail;
413 len -= avail;
414 FInputPos += avail;
415
416 if (FInputPos = FWindowSize) then
417 begin
418 FillFirstHistoryWindow();
419 if (FDoProfile) then
420 GetProfile()
421 else
422 RemoveNoise();
423 RotateHistoryWindows();
424
425 // Rotate halfway for overlap-add
426 //for(i = 0; i < mWindowSize / 2; i++) {
427 for i := 0 to FWindowSize div 2 -1 do
428 FInWaveBuffer[i] := FInWaveBuffer[i + FWindowSize div 2];
429
430 FInputPos := FWindowSize div 2;
431 end;
432 end;
433end;
434
435procedure TNoiseRemoval.FillFirstHistoryWindow;
436var
437 i: Integer;
438begin
439 for i := 0 to FWindowSize-1 do
440 FFFTBuffer[i] := FInWaveBuffer[i];
441 FFFT^.RealFFTf(FFFTBuffer);
442 //for(i = 1; i < (mSpectrumSize-1); i++) {
443 for i := 1 to FSpectrumSize-2 do
444 begin
445 FRealFFTs[0][i] := FFFTBuffer[FFFT^.BitReversed[i] ];
446 FImagFFTs[0][i] := FFFTBuffer[FFFT^.BitReversed[i]+1];
447 FSpectrums[0][i] := FRealFFTs[0][i]*FRealFFTs[0][i] + FImagFFTs[0][i]*FImagFFTs[0][i];
448 FGains[0][i] := FNoiseAttenFactor;
449 end;
450
451 // DC and Fs/2 bins need to be handled specially
452 FSpectrums[0][0] := FFFTBuffer[0]*FFFTBuffer[0];
453 FSpectrums[0][FSpectrumSize-1] := FFFTBuffer[1]*FFFTBuffer[1];
454 FGains[0][0] := FNoiseAttenFactor;
455 FGains[0][FSpectrumSize-1] := FNoiseAttenFactor;
456end;
457
458function Max(A,B: Integer): Integer; inline;
459begin
460 if A>B then
461 Exit(A);
462 Result := B;
463end;
464
465procedure TNoiseRemoval.ApplyFreqSmoothing(ASpec: PSingle);
466var
467 tmp: PSingle;
468 i, j, j0, j1: Integer;
469begin
470 tmp := NewFloatArray(FSpectrumSize);
471 for i := 0 to FSpectrumSize-1 do
472 begin
473 j0 := Max(0, i - FFreqSmoothingBins);
474 j1 := Min(FSpectrumSize-1, i + FFreqSmoothingBins);
475 tmp[i] := 0.0;
476 //for(j = j0; j <= j1; j++)
477 for j := j0 to j1-1 do
478 begin
479 tmp[i] += Aspec[j];
480 end;
481 tmp[i] := tmp[i] / (j1 - j0 + 1);
482 end;
483
484 //for(i = 0; i < mSpectrumSize; i++)
485 for i := 0 to FSpectrumSize-1 do
486 Aspec[i] := tmp[i];
487
488 Freemem(Tmp);
489end;
490
491procedure TNoiseRemoval.GetProfile;
492var
493 start,
494 finish,
495 i, j: Integer;
496 min: Single;
497begin
498 // The noise threshold for each frequency is the maximum
499 // level achieved at that frequency for a minimum of
500 // mMinSignalBlocks blocks in a row - the max of a min.
501
502 start := FHistoryLen - FMinSignalBlocks;
503 finish := FHistoryLen;
504
505
506 for j := 0 to FSpectrumSize-1 do
507 begin
508 min := FSpectrums[start][j];
509 for i := start+1 to finish-1 do
510 if (FSpectrums[i][j] < min) then
511 min := FSpectrums[i][j];
512
513 if (min > FNoiseThreshold[j]) then
514 FNoiseThreshold[j] := min;
515 end;
516
517 FOutSampleCount += FWindowSize div 2; // what is this for? Not used when we are getting the profile?
518end;
519
520procedure TNoiseRemoval.RemoveNoise;
521var
522 center: Integer;
523 start,
524 finish,
525 i,j : Integer;
526 min: Single;
527 out_: Integer;
528begin
529 center := FHistoryLen div 2;
530 start := center - FMinSignalBlocks div 2;
531 finish := start + FMinSignalBlocks;
532
533 // Raise the gain for elements in the center of the sliding history
534 for j := 0 to FSpectrumSize-1 do
535 begin
536 min := FSpectrums[start][j];
537 //for (i = start+1; i < finish; i++) {
538 for i := start+1 to finish-1 do
539 begin
540 if (FSpectrums[i][j] < min) then
541 min := FSpectrums[i][j];
542 end;
543 if (min > FSensitivityFactor * FNoiseThreshold[j]) and (FGains[center][j] < 1.0) then
544 begin
545 if (FbLeaveNoise) then
546 FGains[center][j] := 0.0
547 else
548 FGains[center][j] := 1.0;
549 end
550 else
551 begin
552 if (FbLeaveNoise) then
553 FGains[center][j] := 1.0;
554 end;
555 end;
556
557 // Decay the gain in both directions;
558 // note that mOneBlockAttackDecay is less than 1.0
559 // of linear attenuation per block
560 for j := 0 to FSpectrumSize-1 do
561 begin
562 for i := center+1 to FHistoryLen-1 do
563 begin
564 if (FGains[i][j] < FGains[i - 1][j] * FOneBlockAttackDecay) then
565 FGains[i][j] := FGains[i - 1][j] * FOneBlockAttackDecay;
566 if (FGains[i][j] < FNoiseAttenFactor) then
567 FGains[i][j] := FNoiseAttenFactor;
568 end;
569 for i := center-1 downto 0 do
570 begin
571 if (FGains[i][j] < FGains[i + 1][j] * FOneBlockAttackDecay) then
572 FGains[i][j] := FGains[i + 1][j] * FOneBlockAttackDecay;
573 if (FGains[i][j] < FNoiseAttenFactor) then
574 FGains[i][j] := FNoiseAttenFactor;
575 end;
576 end;
577
578
579 // Apply frequency smoothing to output gain
580 out_ := FHistoryLen - 1; // end of the queue
581
582 ApplyFreqSmoothing(FGains[out_]);
583
584 // Apply gain to FFT
585 //for (j = 0; j < (mSpectrumSize-1); j++) {
586 for j := 0 to FSpectrumSize-2 do
587 begin
588 FFFTBuffer[j*2 ] := FRealFFTs[out_][j] * FGains[out_][j];
589 FFFTBuffer[j*2+1] := FImagFFTs[out_][j] * FGains[out_][j];
590 end;
591 // The Fs/2 component is stored as the imaginary part of the DC component
592 FFFTBuffer[1] := FRealFFTs[out_][FSpectrumSize-1] * FGains[out_][FSpectrumSize-1];
593
594 // Invert the FFT into the output buffer
595 FFFT^.InverseRealFFTf(FFFTBuffer);
596
597 // Overlap-add
598 for j := 0 to FSpectrumSize-2 do
599 begin
600 FOutOverlapBuffer[j*2 ] += FFFTBuffer[FFFT^.BitReversed[j] ] * FWindow[j*2 ];
601 FOutOverlapBuffer[j*2+1] += FFFTBuffer[FFFT^.BitReversed[j]+1] * FWindow[j*2+1];
602 end;
603
604 // Output the first half of the overlap buffer, they're done -
605 // and then shift the next half over.
606
607 if (FOutSampleCount >= 0) then // ...but not if it's the first half-window
608 begin
609 //FOutputTrack->Append((samplePtr)mOutOverlapBuffer, floatSample, mWindowSize / 2);
610 FWriteProc(Self, FOutOverlapBuffer, FWindowSize div 2);
611 end;
612
613
614 FOutSampleCount += FWindowSize div 2;
615 //for(j = 0; j < mWindowSize / 2; j++)
616 for j := 0 to FWindowSize div 2 -1 do
617 begin
618 FOutOverlapBuffer[j] := FOutOverlapBuffer[j + (FWindowSize div 2)];
619 FOutOverlapBuffer[j + (FWindowSize div 2)] := 0.0;
620 end
621end;
622
623procedure TNoiseRemoval.RotateHistoryWindows;
624var
625 last: Integer;
626 i: Integer;
627 lastSpectrum: PSingle;
628 lastGain: PSingle;
629 lastRealFFT: PSingle;
630 lastImagFFT: PSingle;
631begin
632 last := FHistoryLen - 1;
633
634 // Remember the last window so we can reuse it
635 lastSpectrum := FSpectrums[last];
636 lastGain := FGains[last];
637 lastRealFFT := FRealFFTs[last];
638 lastImagFFT := FImagFFTs[last];
639
640 // Rotate each window forward
641 //for(i = last; i >= 1; i--) {
642 for i := last downto 1 do
643 begin
644 FSpectrums[i] := FSpectrums[i-1];
645 FGains[i] := FGains[i-1];
646 FRealFFTs[i] := FRealFFTs[i-1];
647 FImagFFTs[i] := FImagFFTs[i-1];
648 end;
649
650 // Reuse the last buffers as the new first window
651 FSpectrums[0] := lastSpectrum;
652 FGains[0] := lastGain;
653 FRealFFTs[0] := lastRealFFT;
654 FImagFFTs[0] := lastImagFFT;
655end;
656
657procedure TNoiseRemoval.FinishTrack;
658var
659 empty: PSingle;
660 i: Integer;
661begin
662 // Keep flushing empty input buffers through the history
663 // windows until we've output exactly as many samples as
664 // were input.
665 // Well, not exactly, but not more than mWindowSize/2 extra samples at the end.
666 // We'll delete them later in ProcessOne.
667 empty := NewFloatArray(FWindowSize div 2);
668 //for(i = 0; i < mWindowSize / 2; i++)
669 for i := 0 to FWindowSize div 2 -1 do
670 empty[i] := 0.0;
671
672 while (FOutSampleCount < FInSampleCount) do
673 ProcessSamples(FWindowSize div 2, empty);
674
675 Freemem(empty);
676end;
677
678procedure TNoiseRemoval.Cleanup;
679var
680 i: Integer;
681begin
682 FFFT^.EndFFT;
683
684 if (FDoProfile) then
685 ApplyFreqSmoothing(@FNoiseThreshold[0]);
686
687
688 for i := 0 to FHistoryLen-1 do
689 begin
690 FreeMem(FSpectrums[i]);
691 FreeMem(FGains[i]);
692 FreeMem(FRealFFTs[i]);
693 FreeMem(FImagFFTs[i]);
694 end;
695 SetLength(FSpectrums,0);
696 SetLength(FGains,0);
697 SetLength(FRealFFTs,0);
698 SetLength(FImagFFTs,0);
699
700 FreeMem(FFFTBuffer);
701 FreeMem(FInWaveBuffer);
702 FreeMem(FWindow);
703 FreeMem(FOutOverlapBuffer);
704
705 FInited := False;
706end;
707
708function TNoiseRemoval.GetNoiseProfile: TSingleArray;
709begin
710 SetLength(Result, FSpectrumSize);
711 Move(FNoiseThreshold[0], Result[0], FSpectrumSize);
712end;
713
714procedure TNoiseRemoval.SetAttackDecayTime(AValue: Double);
715begin
716 if FAttackDecayTime=AValue then Exit;
717 if AValue < 0.0 then AValue := 0;
718 if AValue > 1.0 then AValue := 1.0;
719 FAttackDecayTime:=AValue;
720end;
721
722procedure TNoiseRemoval.SetFreqSmoothingHz(AValue: Double);
723begin
724 if FFreqSmoothingHz=AValue then Exit;
725 if AValue<0 then AValue:=0;
726 if AValue>1000 then AValue := 1000;
727 FFreqSmoothingHz:=AValue;
728end;
729
730procedure TNoiseRemoval.SetGain(AValue: Double);
731begin
732 if FNoiseGain=AValue then Exit;
733 if AValue > 0 then AValue:=0;
734 if AValue < -48 then AValue := -48;
735
736 FNoiseGain:=AValue;
737end;
738
739procedure TNoiseRemoval.SetNoiseProfile(AValue: TSingleArray);
740begin
741 SetLength(FNoiseThreshold, FSpectrumSize);
742 Move(AValue[0], FNoiseThreshold[0], FSpectrumSize);
743 FHasProfile:=True;
744
745 FDoProfile:=False;
746 Cleanup; // set after FDoProfile so the profile is not processed.
747
748end;
749
750procedure TNoiseRemoval.SetSensitivity(AValue: Double);
751begin
752 if FSensitivity=AValue then Exit;
753 if AValue < -20.0 then AValue:=-20.0;
754 if AValue > 20.0 then AValue := 20.0;
755 FSensitivity:=AValue;
756end;
757
758constructor TNoiseRemoval.Create;
759begin
760 FWindowSize:=2048;
761 FSpectrumSize:= 1 + FWindowSize div 2;
762
763 // loaded prefs
764 FSensitivity := 0.0;
765 FNoiseGain := -24.0;
766 FFreqSmoothingHz := 150.0;
767 FAttackDecayTime:= 0.15;
768 FbLeaveNoise:=False;
769
770 FMinSignalTime := 0.05;
771 FHasProfile := False;
772 FDoProfile := True;
773
774 SetLength(FNoiseThreshold, FSpectrumSize);
775
776end;
777
778destructor TNoiseRemoval.Destroy;
779begin
780 SetLength(FNoiseThreshold, 0);
781 inherited Destroy;
782end;
783
784function TNoiseRemoval.Init(ASampleRate: Integer): Boolean;
785begin
786 FSampleRate:=ASampleRate;
787 Initialize;
788 FInited:=True;
789 Result := True;
790 Reset;
791end;
792
793function TNoiseRemoval.Process(AData: PSingle; ASampleCount: Integer;
794 AGetNoiseProfile: Boolean; AMoreComing: Boolean = False): Boolean;
795begin
796 if not FInited then
797 Raise Exception.Create('TNoiseRemoval is not Inited');
798
799 if not AGetNoiseProfile and not FHasProfile then
800 raise Exception.Create('Tried to remove noise without profile.');
801
802 FDoProfile:=AGetNoiseProfile;
803
804 if FDoProfile and (FTotalRead = 0) then
805 begin
806 Initialize;
807 reset;
808 end;
809
810 Inc(FTotalRead, ASampleCount);
811 ProcessSamples(ASampleCount, AData);
812 Result := True;
813
814 if AMoreComing then
815 Exit;
816
817 if {AGetNoiseProfile or }FDoProfile then
818 begin
819 FHasProfile:=True;
820 Cleanup; // triggers the data in FNoiseThreshold to be processed
821
822 // must be set after Cleanup() is called
823 //FDoProfile:=False;
824
825 //Initialize;
826 FTotalRead:=0;
827 end;
828
829 FHasProfile:=True;
830 FDoProfile := False;
831
832end;
833
834procedure TNoiseRemoval.Flush;
835begin
836 if not FInited then
837 Exit;
838
839 FinishTrack;
840 Cleanup; // Sets FInited to False
841end;
842
843{ TFFT }
844
845function TFFT.InitializeFFT(FFTLen: Integer): PFFT;
846var
847 i: Integer;
848 temp: Integer;
849 mask: Integer;
850begin
851 Result := New(PFFT);
852 if Result = nil then
853 raise EOutOfMemory.Create('Error allocating memory for FFT');
854
855
856 with Result^ do begin
857
858 {*
859 * FFT size is only half the number of data points
860 * The full FFT output can be reconstructed from this FFT's output.
861 * (This optimization can be made since the data is real.)
862 *}
863 Points := FFTLen div 2;
864
865 SetLength(FPCSinTable, 2*Points);
866 SinTable:=@FPCSinTable[0];
867
868 BitReversed := Getmemory(Points*SizeOf(Integer));
869 if BitReversed = nil then
870 raise EOutOfMemory.Create('Error allocating memory for BitReversed.');
871
872 for i := 0 to Points-1 do
873 begin
874 temp:=0;
875 mask := Points div 2;
876 while mask > 0 do
877 begin
878 //for(mask=h->Points/2;mask>0;mask >>= 1)
879 // temp=(temp >> 1) + (i&mask ? h->Points : 0);
880 temp := (temp shr 1);
881 if (i and mask) <> 0 then
882 temp := temp + Points;
883 //else temp := temp + 0; // why would you do that?
884 mask := mask shr 1;
885 end;
886
887 BitReversed[i]:=temp;
888 end;
889
890 for i := 0 to Points-1 do
891 begin
892 SinTable[BitReversed[i] ]:= -sin(2*PI*i/(2*Points));
893 SinTable[BitReversed[i]+1]:= -cos(2*PI*i/(2*Points));
894 end;
895
896{$ifdef EXPERIMENTAL_EQ_SSE_THREADED}
897 // new SSE FFT routines work on live data
898 for(i=0;i<32;i++)
899 if((1<<i)&fftlen)
900 h->pow2Bits=i;
901{$endif}
902
903 end; // with Result^
904
905
906end;
907
908procedure TFFT.EndFFT;
909begin
910 if Points>0 then
911 begin
912 Freemem(BitReversed);
913 SetLength(FPCSinTable, 0);
914 SinTable:=nil;
915 end;
916
917 Dispose(PFFT(@Self));
918end;
919
920function TFFT.GetFFT(FFTLen: Integer): PFFT;
921var
922 h: Integer = 0;
923 n: Integer;
924begin
925 n := fftlen div 2;
926
927 while (h<MAX_HFFT) and (FFTArray[h] <> nil) and (n <> FFTArray[h]^.Points) do
928 begin
929 if (h<MAX_HFFT) then
930 begin
931 if(FFTArray[h] = nil) then
932 begin
933 FFTArray[h] := InitializeFFT(fftlen);
934 FFTLockCount[h] := 0;
935 end;
936 Inc(FFTLockCount[h]);
937 Exit(FFTArray[h]);
938 end
939 else begin
940 // All buffers used, so fall back to allocating a new set of tables
941 Exit(InitializeFFT(fftlen));
942 end;
943 Inc(h);
944 end;
945end;
946
947procedure TFFT.ReleaseFFT;
948var
949 h: Integer = 0;
950begin
951
952 while (h<MAX_HFFT) and (FFTArray[h] <> @Self) do
953 begin
954 if(h<MAX_HFFT) then
955 begin
956 Dec(FFTLockCount[h]);
957 end
958 else
959 begin
960 EndFFT;
961 end;
962 Inc(h);
963 end;
964end;
965
966procedure TFFT.InverseRealFFTf(buffer: PSingle);
967var
968 A, B: PSingle;
969 sptr: PSingle;
970 endptr1, endptr2: PSingle;
971 br1: PInteger;
972 HRplus,HRminus,HIplus,HIminus: Single;
973 v1,v2,sin,cos: Single;
974 ButterfliesPerGroup: Integer;
975begin
976
977 ButterfliesPerGroup:=Points div 2;
978
979 //* Massage input to get the input for a real output sequence. */
980 A:=@buffer[2];
981 B:=@buffer[Points*2-2];
982 br1:=@BitReversed[1];
983 while(A<B) do
984 begin
985 sin:=SinTable[br1^];
986 cos:=SinTable[br1[1]];
987 //HRplus = (HRminus = *A - *B ) + (*B * 2);
988 HRminus:=A^-B^;
989 HRplus:=HRminus+ (B^ *2);
990
991 //HIplus = (HIminus = *(A+1) - *(B+1)) + (*(B+1) * 2);
992 HIminus:=A[1]-B[1];
993 HIplus:=HIminus+(B[1] *2);
994
995 v1 := (sin*HRminus + cos*HIplus);
996 v2 := (cos*HRminus - sin*HIplus);
997 A^ := (HRplus + v1) * single(0.5);
998 B^ := A^ - v1;
999 A[1] := (HIminus - v2) * single(0.5);
1000 B[1] := A[1] - HIminus;
1001
1002 A+=2;
1003 B-=2;
1004 Inc(br1);
1005 end;
1006 //* Handle center bin (just need conjugate) */
1007 A[1] :=-A[1];
1008 {* Handle DC bin separately - this ignores any Fs/2 component
1009 buffer[1]=buffer[0]=buffer[0]/2;*}
1010 //* Handle DC and Fs/2 bins specially */
1011 //* The DC bin is passed in as the real part of the DC complex value */
1012 //* The Fs/2 bin is passed in as the imaginary part of the DC complex value */
1013 //* (v1+v2) = buffer[0] == the DC component */
1014 //* (v1-v2) = buffer[1] == the Fs/2 component */
1015 v1:=0.5*(buffer[0]+buffer[1]);
1016 v2:=0.5*(buffer[0]-buffer[1]);
1017 buffer[0]:=v1;
1018 buffer[1]:=v2;
1019
1020 {*
1021 * Butterfly:
1022 * Ain-----Aout
1023 * \ /
1024 * / \
1025 * Bin-----Bout
1026 *}
1027
1028 endptr1:=@buffer[Points*2];
1029
1030 while(ButterfliesPerGroup>0) do
1031 begin
1032 A:=buffer;
1033 B:=@buffer[ButterfliesPerGroup*2];
1034 sptr:=@SinTable[0];
1035
1036 while(A<endptr1) do
1037 begin
1038 sin:=sptr^; Inc(sptr); // *(sptr++);
1039 cos:=sptr^; Inc(sptr); // *(sptr++);
1040 endptr2:=B;
1041 while(A<endptr2) do
1042 begin
1043 v1:=B^*cos - B[1]*sin;
1044 v2:=B^*sin + B[1]*cos;
1045 B^ := (A^+v1)*Single(0.5);
1046 A^ := B^ - v1; Inc(A); Inc(B); //*(A++)=*(B++)-v1;
1047 B^ := (A^ + v2)* Single(0.5); //*B=(*A+v2)*(fft_type)0.5;
1048 A^ := B^ - v2; Inc(A); Inc(B); //*(A++)=*(B++)-v2;
1049 end;
1050 A:=B;
1051 B := @B[ButterfliesPerGroup*2];
1052 end;
1053 ButterfliesPerGroup := ButterfliesPerGroup shr 1;
1054 end;
1055end;
1056
1057procedure TFFT.CleanupFFT;
1058var
1059 h: Integer;
1060begin
1061
1062 for h :=0 to MAX_HFFT-1do begin
1063 if((FFTLockCount[h] <= 0) and (FFTArray[h] <> nil)) then
1064 begin
1065 FFTArray[h]^.EndFFT;
1066 FFTArray[h] := nil;
1067 end;
1068 end;
1069end;
1070
1071procedure TFFT.RealFFTf(buffer: PSingle);
1072var
1073 A, B: PSingle;
1074 sptr: PSingle;
1075 endptr1, endptr2: PSingle;
1076 br1, br2: PInteger;
1077 HRplus,HRminus,HIplus,HIminus: Single;
1078 v1,v2,sin_,cos_: Single;
1079 ButterfliesPerGroup: Integer;
1080begin
1081 ButterfliesPerGroup:=Points div 2;
1082
1083 {*
1084 * Butterfly:
1085 * Ain-----Aout
1086 * \ /
1087 * / \
1088 * Bin-----Bout
1089 *}
1090
1091 endptr1:=buffer+Points*2;
1092
1093 while(ButterfliesPerGroup>0) do
1094 begin
1095 A:=buffer;
1096 B:=buffer+ButterfliesPerGroup*2;
1097 sptr:=@SinTable[0];
1098
1099 while(A<endptr1) do
1100 begin
1101 sin_:=sptr^;
1102 cos_ := sptr[1];
1103 endptr2:=B;
1104 while(A<endptr2) do
1105 begin
1106 v1 := B^ * cos_ + B[1] * sin_; //v1=*B*cos + *(B+1)*sin;
1107 v2 := B^ * sin_ - B[1] * cos_; //v2=*B*sin - *(B+1)*cos;
1108 B^ := A^+v1; //*B=(*A+v1);
1109 A^ := B^-2*v1; Inc(A); Inc(B); //*(A++)=*(B++)-2*v1;
1110
1111 B^ := A^-v2; //*B=(*A-v2);
1112 A^ := B^+2*v2; Inc(A); Inc(B); //*(A++)=*(B++)+2*v2;
1113 end;
1114 A:=B;
1115 B:=B+ButterfliesPerGroup*2;
1116 sptr:=sptr+2;
1117 end;
1118
1119 ButterfliesPerGroup := ButterfliesPerGroup shr 1;
1120 end;
1121
1122 //* Massage output to get the output for a real input sequence. */
1123 br1:=@BitReversed[1]; // is this wrong? Should be @BitReversed[0] ; ?
1124 br2:=@BitReversed[Points-1];
1125
1126 while(br1<br2) do
1127 begin
1128 sin_:=SinTable[br1[0]];
1129 cos_:=SinTable[br1[1]];
1130 A:=@buffer[br1^];
1131 B:=@buffer[br2^];
1132 //HRplus = (HRminus = *A - *B ) + (*B * 2);
1133 HRminus := A^ - B^;
1134 HRplus := HRminus + (B^ * 2);
1135
1136 //HIplus = (HIminus = *(A+1) - *(B+1)) + (*(B+1) * 2);
1137 HIminus := A[1] - B[1];
1138 HIplus := HIminus + (B[1] * 2);
1139
1140 v1 := (sin_*HRminus - cos_*HIplus);
1141 v2 := (cos_*HRminus + sin_*HIplus);
1142 A^ := (HRplus + v1) * single(0.5);
1143 B^ := A^ - v1;
1144 A[1] := (HIminus + v2) * single(0.5);
1145 B[1] := A[1] - HIminus;
1146
1147 Inc(br1);
1148 Dec(br2);
1149 end;
1150 //* Handle the center bin (just need a conjugate) */
1151 A:=buffer+br1[1];
1152 A^ := -A^;
1153 {* Handle DC bin separately - and ignore the Fs/2 bin
1154 buffer[0]+=buffer[1];
1155 buffer[1]=(fft_type)0;*}
1156 ///* Handle DC and Fs/2 bins separately */
1157 ///* Put the Fs/2 value into the imaginary part of the DC bin */
1158 v1:=buffer[0]-buffer[1];
1159 buffer[0]+=buffer[1];
1160 buffer[1]:=v1;
1161
1162end;
1163
1164procedure TFFT.ReorderToTime(Buffer: PSingle; TimeOut: PSingle);
1165var
1166 i: Integer;
1167begin
1168 // Copy the data into the real outputs
1169 //for(int i=0;i<hFFT->Points;i++) {
1170 for i := 0 to Points-1 do
1171 begin
1172 TimeOut[i*2 ]:=buffer[BitReversed[i] ];
1173 TimeOut[i*2+1]:=buffer[BitReversed[i]+1];
1174 end;
1175end;
1176
1177procedure TFFT.ReorderToFreq(Buffer: PSingle; RealOut: PSingle; ImagOut: PSingle
1178 );
1179var
1180 i: Integer;
1181begin
1182 // Copy the data into the real and imaginary outputs
1183 //for(int i=1;i<hFFT->Points;i++)
1184 for i := 1 to Points-1 do
1185 begin
1186
1187 RealOut[i]:=buffer[BitReversed[i] ];
1188 ImagOut[i]:=buffer[BitReversed[i]+1];
1189 end;
1190 RealOut[0] := buffer[0]; // DC component
1191 ImagOut[0] := 0;
1192 RealOut[Points] := buffer[1]; // Fs/2 component
1193 ImagOut[Points] := 0;
1194end;
1195
1196initialization
1197 FillChar(FFTArray, SizeOf(FFTArray), 0);
1198 FillChar(FFTLockCount, SizeOf(FFTLockCount), 0);
1199
1200end.
1201
Note: See TracBrowser for help on using the repository browser.