Tweet
#ifndef PI
#define PI 3.14159265358979323846264338327950288419716939937510
#endif
#ifndef PI2
#define PI2 (2.0*PI)
#endif
//-----------------------------------------------------------------------------
// Data Types
//
typedef struct // Hilbert Filter Object Definition
{
int N; // Total Filter Size (N), must be odd (e.g. 15, 31, 63, ..)
int pos; // Current Filter Position
int nC; // Total Number of Stored Coefficients (0.0 and symmetry excluded)
REAL *buffer; // Double-Buffer to avoid modulo indexing (linearized indexing)
REAL *coeffs; // Hilbert Coefficients (odd, positive)
} HilbertFilter;
//-----------------------------------------------------------------------------
// Functions
//
//-----------------------------------------------------------------------------
// Initialize Hilbert Transformation and Hilbert Filter Object generation
// [IN] int N: Filter size (must be odd size)
// [IN] DWORD dwSR: Sample rate in Hz, if 0 -> apply no band pass filter (=allpass)
// [IN] REAL f_min: Lower corner frequency for band pass filter in Hz
// [IN] REAL f_max: Upper corner frequency for band pass filter in Hz
// [IN] BOOL bHamming: 1=apply Hamming-Windowing, 0=no windowing
//
// [OUT] HilbertFilter*: Pointer to Hilbert-Object or NULL if not enough memory available
HilbertFilter* hilbert_init(int N, DWORD dwSR, REAL f_min, REAL f_max, BOOL bHamming)
{
HilbertFilter *hf;
int i, M, n;
REAL hamming=1.0;
REAL bp=1.0;
REAL hic;
REAL bw_low;
REAL bw_high;
if (N % 2 == 0) // Make N odd if required
N++;
hf = MALLOC(sizeof(HilbertFilter));
ASSERT(hf);
hf->N = N;
hf->nC = (N - 1) / 4 + ((N - 1) / 2 % 2 != 0 ? 1 : 0); // Nur ungerade Indizes
hf->coeffs = MALLOC(hf->nC * sizeof(REAL));
ASSERT(hf->coeffs);
hf->buffer = CALLOC(2 * N, sizeof(REAL));
ASSERT(hf->buffer);
if (dwSR != 0)
{
bw_low = PI2 * f_min / (REAL)dwSR;
bw_high = PI2 * f_max / (REAL)dwSR;
}
hf->pos = 0;
M = N / 2;
n = 0;
for (i = 1; i <= M; i += 2) // for all odd index i
{
// Ideal Hilbert Coefficients for index i (allpass)
hic = 2.0 / (PI * i);
// Hamming-Windowing Coefficients (optional)
if (bHamming)
hamming = 0.54 + 0.46 * cos((PI2 * i) / (N - 1));
// Band Pass Filter Coefficients (optional)
if (dwSR != 0)
bp = cos(i * bw_low) - cos(i * bw_high);
// Combined Hilbert Coefficients (all-in-one)
hf->coeffs[n++] = hic * bp * hamming;
} // for i
hf->nC = n;
return hf;
}
//-----------------------------------------------------------------------------
// Perform a Hilbert Transformation and get I/Q channel results for the input sample
// [IN] HilbertFilter *hf: Pointer to Hilbert Object
// [IN] REAL input: Input sample (one sample)
// [OUT] REAL *quadrature: Quadrature output (one sample)
//
// [OUT] REAL Inphase output (one Sample)
REAL hilbert_process(HilbertFilter *hf, REAL input, REAL *quadrature)
{
int i, n, N, M;
REAL q, inphase, *window;
ASSERT(hf);
ASSERT(quadrature);
N = hf->N;
M = N / 2;
// Lese-Basis für das älteste Fenster
window = &hf->buffer[hf->pos];
// Double Buffering: Sample an zwei Positionen schreiben
window[0] = input;
window[N] = input;
q = 0.0;
// Symmetrie-Optimierung: Nur ungerade n nutzen
for (i = 0; i < hf->nC; i++)
{
n = 2 * i + 1;
q += hf->coeffs[i] * (window[M + n] - window[M - n]);
} // for i
*quadrature = q;
inphase = window[(N - 1) / 2]; // Delayed center sample
hf->pos = (hf->pos + 1) % N; // Modulo nur für den Pointer, nicht pro Multiplikation
return inphase;
}
//-----------------------------------------------------------------------------
// Hilbert Objekt freigeben
// [IN] HilbertFilter *hf: Zeiger auf Hilbert-Objekt
// [OUT] none
void hilbert_free(HilbertFilter *hf)
{
ASSERT(hf);
ASSERT(hf->buffer);
ASSERT(hf->coeffs);
FREE(hf->buffer);
FREE(hf->coeffs);
FREE(hf);
}
The posted code for unwrapping the phase doesn't work. I have taken another (simple) solution, which works.// Phase Unwrapping (Phase linearising)
// [OUT] REAL *pOutput: Pointer to Output Phase Data (Radiants)
// [IN] REAL *pInput: Pointer to Input Phase Data (Radiants)
// [IN] int n: Number of Phase Data
//
void unwrapPhase1D(REAL *pOutput, REAL *pInput, int n)
{
REAL wrap_count, diff;
int i;
ASSERT(pOutput);
ASSERT(pInput);
if (n <= 0)
return;
wrap_count = 0;
pOutput[0] = pInput[0];
for (i = 1; i < n; i++)
{
// Current difference
diff = pInput[i] - pInput[i-1];
// Detection of phase jumps > pi or < -pi
if (diff > PI) {
wrap_count -= 1;
} else if (diff < -PI) {
wrap_count += 1;
}
pOutput[i] = pInput[i] + wrap_count * PI2;
} // for i
}
Comment