|
ru.algorithms
- RU.ALGORITHMS ----------------------------------------------------------------
From : Alexander Rebrov 2:5000/115.17 11 Nov 2001 23:43:36
To : Aleksey V Vaneev
Subject : Re: FFT
--------------------------------------------------------------------------------
09 ноября 2001 21:37, you wrote to me:
AV>>> она веpхи (imho, выше 10кГц) pубает почем зpя.
AR>> 4-х точечная кубическая? Если да, то она, имхо, ничего
AR>> не обpезает.
AV> Если бы это было так, то зачем было бы использовать формулы более
AV> высокого порядка? :) А вообще, она, imho, дает т.н. "muddy sound".
AV> Чистый, но нечеткий. Hаверное, по фазам перекос дает.
Ты кубическую так? делаешь? Я попробовал результат работы в виде графика
посмотреть. По-моему Hикаких перекосов. А то что звук Hе совсем четкий - так
это точек все же маловато.
===>
fPosition := fPosition + fRatio;
nPosition := Trunc (fPosition);
p := fPosition - nPosition;
q := (1.0 - p);
if nOldPos <> nPosition Then
Begin
pBPos := pBPos + (nPosition - nOldPos);
If pBPos > nInSize - 2 Then
Begin
Result:=samp;
Exit;
End;
ABCD := pIn^[pBPos] - pIn^[pBPos-1] + pIn^[pBPos+1] - pIn^[pBPos+2];
CB := pIn^[pBPos+1] - pIn^[pBPos];
B := pIn^[pBPos];
nOldPos := nPosition;
End;
tmp := B + p*(0.25*q*ABCD+CB);
if tmp > 32767 Then pOut^[samp] := 32767 Else
if tmp < -32767 Then pOut^[samp] := -32767 Else
pOut^[samp] := Round (tmp);
Inc (Samp);
===>
AR>> А вообще есть windowed-sinc с окном Кайзеpа. Качество идеальное.
AR>> Там этих точек хоть миллион задавай. Пpобовал пpимеpно со 170 до
AR>> 44100 интеpполиpовать - абсолютно чистый звук (эс-ноодни низы).
AV> А у тебя есть исходники, может, зашлешь мылом?
А вот sinc иHтерполяция. Тут кто-то тоже Hасчет этого базарит, поэтому сюда
закиHу:
===>
Const
MaxInterpFilterSize = 16384;
var
// Массив значений фильтра (т.н. impulse response)
InterpFilterImp: Array [0..MaxInterpFilterSize-1] Of Single;
// Массив разниц: FilterImpDeltas[i]:=FilterImp[i+1]-FilterImp[i]
// Служит для быстрой линейной интерполяции фильтра(!), т.е.
// к интерполирующему фильтру применяется линейная интерполяция.
InterpFilterImpDeltas: Array [0..MaxInterpFilterSize-1] Of Single;
// Количество точек интерполяции в одну сторону
InterpNumPoints: Integer;
// Количество значений функции в одном периоде
InterpFilterPeriodSize: Integer;
// Количество значений фильтра
InterpFilterSize: Integer;
{ NumPoints - количество точек интерполяции в одну сторону. Так как при
интерполяции
половина значений берется по левую сторону от текущего сэмпла,
а половина - по правую, то общее количество точек интерполяции
равно NumPoints*2. В связи с этим и количество периодов функции
фильтра равно NumPoints.
PeriodSize - количество значений функции фильтра в одном периоде. Чем больше
значений - тем выше точность интерполяции. Общее кол-во значений
массива равно NumPoints*PeriodSize.
FrqRoll - коэффициент ФЧ, должен быть от 0 до 1; определает частоту среза.
ачинает обрезать частоты свыше
FrqRoll*оваяЧастотаДикретизации/2.
Beta - сам не понял, но это параметр окна Кайзера. Вроде как влияет на
соотношение сигнал/шум. Вот описание на англ.
Beta trades the rejection of the lowpass filter against the"
transition width from passband to stopband. Larger Beta means"
a slower transition and greater stopband rejection. See Rabiner"
and Gold (Th. and App. of DSP) under Kaiser windows for more about"
Beta. The following table from Rabiner and Gold (p. 101) gives
some"
feel for the effect of Beta:"
BETA D PB RIP SB RIP"
2.120 1.50 +-0.27 -30"
3.384 2.23 0.0864 -40"
4.538 2.93 0.0274 -50"
5.658 3.62 0.00868 -60"
6.764 4.32 0.00275 -70"
7.865 5.0 0.000868 -80"
8.960 5.7 0.000275 -90"
10.056 6.4 0.000087 -100"
Above, ripples are in dB, and the transition band width is "
approximately D*Fs/N, where Fs = sampling rate, "
and N = window length. PB = 'pass band' and SB = 'stop band'."
Alternatively, D is the transition widths in bins given a"
length N DFT (i.e. a window transform with no zero padding."
}
Function CreateFilter (NumPoints, PeriodSize: Integer; FrqRoll, Beta: Single):
Boolean;
Function Izero (X: Double): Double;
Const
IzeroEPSILON = 1E-21; // Max error acceptable in Izero
Var
Sum, U, HalfX, Temp: Double;
N: Integer;
Begin
Sum := 1; U := 1; N := 1;
HalfX := X * 0.5;
Repeat
Temp := HalfX / N;
Inc (N);
Temp := Temp * Temp;
U := U * Temp;
Sum := Sum + U;
Until U < IzeroEPSILON * Sum;
Result := Sum;
End;
Var
IBeta, Temp1, Temp2, Inm1, alpha, alpha1, alpha2: Double;
I: Integer;
Begin
InterpNumPoints := NumPoints;
InterpFilterPeriodSize := PeriodSize;
InterpFilterSize := NumPoints * PeriodSize;
// Проверяем размер массива
If InterpFilterSize > MaxInterpFilterSize Then
Begin
Result := False;
Exit;
End;
// Calculate ideal lowpass filter impulse response coefficients:
alpha := sin (pi / (InterpFilterSize - 1)) / sin (2 * pi / (InterpFilterSize
- 1));
alpha := alpha * alpha;
alpha1 := (alpha + 1) / 2;
alpha2 := alpha / 2;
InterpFilterImp[0] := FrqRoll;
IBeta := 1.0 / Izero (Beta);
Inm1 := 1.0 / (InterpFilterSize-1);
For I:=1 To InterpFilterSize-1 Do
Begin
// по идее сначала нужно посчитать все значения InterpFilterImp,
// затем найти максимальное, вычислить коэффициент, при умножении на
// который это максимальное будет равно еденице, а затем перемножить
// все значения фильтра на этот коэффициент и затем уж посчитать
// InterpFilterImpDeltas. о на практике оказалось, что оно вроде как
// изначально все правильно.
Temp1 := PI * I / PeriodSize;
Temp2 := I * Inm1;
// Analog sinc-функции, частота среза = FrqRoll
// В принципе, функция sinc(x)=sin(pi*x)/(pi*x) - это та же самая
// интерполяционная функция Лагранжа, только проще.
InterpFilterImp[i] :=
sin (Temp1 * FrqRoll) / Temp1 *
// Calculate and Apply Kaiser window to ideal lowpass filter.
// Note: last window value is IBeta which is NOT zero.
// You're supposed to really truncate the window here, not ramp
// it to zero. This helps reduce the first sidelobe.
Izero (Beta * sqrt (1.0 - temp2 * temp2)) * IBeta;
InterpFilterImpDeltas[i-1]:=InterpFilterImp[i]-InterpFilterImp[i-1];
End;
// Последний коэффициент не интерполируется
InterpFilterImpDeltas[InterpFilterSize-1] :=
-InterpFilterImp[InterpFilterSize-1];
Result := True;
End;
{ ****************************************************** SincInterpolation ***
}
Function SincInterpolation (fRatio: Single; pIn, pOut: PSmallIntArray; nInSize:
Integer; nOutSize: Integer; FLF: Boolean): Integer;
Var
samp, x, tmp, Hoi: Integer;
v, t, Ho, dh: double;
nPosition: Integer; // index of first sample
fPosition, p: Single; // floating point position accumulator
sum1: integer;
sum2: double;
points: integer;
PC: integer;
Begin
points:=8;
// В принципе, создать фильтр нужно только один раз.
CreateFilter (points, 512, 1, 9.5);
// Initial input array element pointer
nPosition:= 0; // index of first sample
fPosition:= 0.0; // floating point position accumulator
// Filter sampling period
dh := InterpFilterPeriodSize;
// for all output samples
pOut^[0]:=pIn^[0];
samp := 0; // =1
p:=0;
If FLF Then
While samp < nOutSize Do // с фильтром Ч
Begin
// THE EQUATION
v:=0;
x:=nPosition; Ho := p * dh;
PC:=0;
While Ho < InterpFilterSize Do
Begin
if x<0 Then Break;
Hoi:=Trunc(Ho);
v:=v+pIn^[x]*(InterpFilterImp[Hoi]+(Ho-Hoi)*InterpFilterImpDeltas[H
oi]);
Ho:=Ho+dh;
x:=x-1;
Inc (PC);
End;
If PC<Points then
PC:=3;
PC:=0;
x:=nPosition+1; Ho := (1-p) * dh;
// If (1-p)=0 Then Ho:=Ho+dh;
While Ho < InterpFilterSize Do
Begin
if x>=nInSize Then Break;
Hoi:=Trunc(Ho);
v:=v+pIn^[x]*(InterpFilterImp[Hoi]+(Ho-Hoi)*InterpFilterImpDeltas[H
oi]);
Ho:=Ho+dh;
x:=x+1;
Inc (PC);
End;
If PC<Points then
PC:=3;
tmp:=Round(v);
// SATURATE AND WRITE RESULTS
if tmp > 32767 Then pOut^[samp] := 32767 Else
if tmp < -32767 Then pOut^[samp] := -32767 Else
pOut^[samp] := Round (tmp);
Inc (Samp);
// increment position by the ratio of # of in samples to # of out
samples
fPosition := fPosition + fRatio;
nPosition := Trunc (fPosition);
p := fPosition - nPosition;
If nPosition > nInSize Then Break;
End
Else
While samp < nOutSize Do
Begin
// THE EQUATION
v:=0;
x:=nPosition; Ho := p * dh;
While Ho < InterpFilterSize Do
Begin
if x<0 Then Break;
Hoi:=Trunc(Ho);
v:=v+pIn^[x]*InterpFilterImp[Hoi];
Ho:=Ho+dh;
x:=x-1;
End;
x:=nPosition+1; Ho := (1-p) * dh;
If (1-p)=0 Then Ho:=Ho+dh;
While Ho < InterpFilterSize-1 Do
Begin
if x>=nInSize Then Break;
Hoi:=Trunc(Ho);
v:=v+pIn^[x]*InterpFilterImp[Hoi];
Ho:=Ho+dh;
x:=x+1;
End;
tmp:=Round(v);
// SATURATE AND WRITE RESULTS
if tmp > 32767 Then pOut^[samp] := 32767 Else
if tmp < -32767 Then pOut^[samp] := -32767 Else
pOut^[samp] := Round (tmp);
Inc (Samp);
// increment position by the ratio of # of in samples to # of out
samples
fPosition := fPosition + fRatio;
nPosition := Trunc (fPosition);
p := fPosition - nPosition;
If nPosition > nInSize Then Break;
End;
Result:=Samp;
End;
===>
Alexander
---
* Origin: (2:5000/115.17)
Вернуться к списку тем, сортированных по: возрастание даты уменьшение даты тема автор
| Тема: |
Автор: |
Дата: |
|
FFT |
Aleksey V. Vaneev |
03 Nov 2001 00:06:10 |
 Re: FFT |
Alexander Rebrov |
03 Nov 2001 22:14:35 |
  FFT |
Aleksey V. Vaneev |
04 Nov 2001 11:02:05 |
 FFT |
Alex Astafiev |
05 Nov 2001 15:58:26 |
|
FFT |
Aleksey V. Vaneev |
06 Nov 2001 09:23:00 |
 Re FFT |
Alexander Rebrov |
07 Nov 2001 21:12:49 |
|
Re FFT |
Aleksey V. Vaneev |
09 Nov 2001 22:37:01 |
|
Re: FFT |
Alexander Rebrov |
11 Nov 2001 23:43:36 |
|
FFT |
Aleksey V. Vaneev |
12 Nov 2001 10:37:59 |
|
Re: FFT |
Alexander Rebrov |
06 Nov 2001 12:36:07 |
|
FFT |
Alex Astafiev |
08 Nov 2001 14:49:10 |
|
FFT |
Daniel Kamperov |
09 Nov 2001 16:56:32 |
|
FFT |
Ђ¤аҐ© |
04 Nov 2001 22:18:59 |
|
FFT |
Evgenij Masherov |
05 Nov 2001 11:20:02 |
|
FFT |
Aleksey V. Vaneev |
06 Nov 2001 09:17:48 |
|
FFT |
Evgenij Masherov |
06 Nov 2001 11:59:59 |
|
FFT |
Alex Astafiev |
06 Nov 2001 09:35:00 |
|
FFT |
…ўЈҐЁ© Њ иҐа®ў |
05 Nov 2001 22:30:36 |
|
FFT |
Aleksey V. Vaneev |
08 Nov 2001 17:27:00 |
|
FFT |
Evgenij Masherov |
09 Nov 2001 10:48:01 |
|
FFT |
Aleksey V. Vaneev |
09 Nov 2001 22:49:40 |
|
FFT |
Evgenij Masherov |
12 Nov 2001 11:03:58 |
|
|
