Я пытаюсь написать преобразование Гильберта с нуля, но не использую какие-либо встроенные библиотеки, кроме fft и ifft.Я не математик по профессии, но я нашел эти два алгоритма онлайн для преобразования Гильберта, один в C и один в MATLAB.Я пытался реализовать их оба, но ни один из них не дал мне того же результата, что и Гильберт SciPy.Я наверняка сделал что-то не так в своей реализации.Любое понимание будет с благодарностью.
Первая реализация: (с веб-сайта MATLAB) Гильберт использует четырехэтапный алгоритм:
Рассчитать БПФ входной последовательности, сохраняя результат в видеvector x.
Создать вектор h, элементы которого h(i) имеют значения:
1 для i = 1, (n/2)+1
2 для i = 2, 3, ... , (n/2)
0 для i = (n/2)+2, ... , n
Рассчитать поэлементное произведение x и h.
Рассчитать обратное БПФ последовательности, полученной на шаге 3, и вернуть первое n элементов результата.
Моя попытка:
def generate_array(n):
a = np.hstack((np.full(n//2+1, 2), np.zeros(n//2-1)))
a[[0, n//2]] = 1
return a
def hilbert_from_scratch_2(u):
fft_result = fft(u) #scipy fft
n = len(u)
to_multiply = generate_array(n)
result = np.multiply(n,to_multiply)
return ifft(result) #scipy ifft
Вторая реализация: (https://www.cfa.harvard.edu/~spaine/am/download/src/transform.c)
void hilbert(double *z, unsigned long n)
{
double x;
unsigned long i, n2;
n2 = n << 1;
/*
* Compute the (bit-reversed) Fourier transform of z.
*/
fft_dif(z, n);
/*
* Form the transform of the analytic sequence by zeroing
* the transform for negative time, except for the (N/2)th.
* element. Since z is now in bit-reversed order, this means
* zeroing every other complex element. The array indices of
* the elements to be zeroed are 6,7,10,11...etc. (The real
* and imaginary parts of the (N/2)th element are in z[2] and
* z[3], respectively.)
*/
for (i = 6; i < n2; i += 4) {
z[i] = 0.;
z[i+1] = 0.;
}
/*
* The 0th and (N/2)th elements get multiplied by 0.5. Test
* for the trivial 1-point transform, just in case.
*/
z[0] *= 0.5;
z[1] *= 0.5;
if (n > 1) {
z[2] *= 0.5;
z[3] *= 0.5;
}
/*
* Compute the inverse transform.
*/
ifft_dit(z, n);
/*
* Normalize the array. The factor of 2 is left over from
* forming the transform in the time domain.
*/
x = 2. / (double)n;
for (i = 0; i < n2; ++i)
z[i] *= x;
return;
} /* hilbert() */
Моя попытка:
def hilbert_from_scratch(signal):
fast_ft = fft(signal) #scipy fft
for i in range(6,len(signal),4):
fast_ft[i] = 0
fast_ft[i+1] = 0
fast_ft[0] = fast_ft[0]*.5
fast_ft[1] = fast_ft[1]*.5
if(len(fast_ft) > 1):
fast_ft[2] = fast_ft[2]*.5
fast_ft[3] = fast_ft[3]*.5
inverse_fft = ifft(fast_ft) #scipy ifft
x = 2 / len(signal)
for i in range(0,len(signal),1):
inverse_fft[i] = inverse_fft[i]*x
return inverse_fft
Любая идея о том, почему ни один из них не дает такой же результат, как у hilbert SciPy, будет принята с благодарностью.