Néha fent, néha lent – avagy hogyan határozzuk meg a pillanatnyi frekvenciát? – Ups and downs, or how to determine the instantaneous frequency?

A kutatás nemcsak fáradtságos munka után elért elvárt erdemények sorozatából áll – amint ez nemrég kiderült.

Az általam leadott kutatási tervben szerepelt a pillanatnyi frekvencia kiszámítási módjainak a vizsgálata. Az első negyedév elején neki is fogtam ennek a feladatnak, majd a Hilbert-transzformáció és az empirikus AM-FM felbontás tanulmányozása után elkezdtem gondolkodni egy matematikai modell megalkotásán, amely hozzásegítsen a pillanatnyi frekvencia alternatív módon való kiszámításához. Ennek érdekében a módusfüggvényeket, amelyek pillanatnyi frekvenciáját meghatározni akartam, a következő formájúaknak tekintettem ([1] alapján):  \displaystyle c(t)=a(t)\sin(2\pi\omega(t)). Mivel ez egy szinuszoidális jelforma, joggal feltételezhetjük, hogy létezik egy “kiegészítő” párja, amely koszinuszoidális (jelölje ezt  \displaystyle c_q(t)), ketten pedig egy változó amplitúdójú és körfrekvenciájú körmozgás két külön tengelyre való levetítései. A koszinuszoidális jel a szinuszoidális jel kvadratúrája, vagyis  \displaystyle \pi/2-vel késik. Idáig az ötlet megegyezik a Hilbert-transzformáció és az empirikus AM-FM felbontás alapötletével: kiszámolni valamilyen módon a módusfüggvény kvadratúráját, majd a két jel segítségével meghatározni a pillanatnyi fázist, majd a pillanatnyi frekvenciát.

A mostani bejegyzésben az egyik zsákutcát szeretném bemutatni, ami a pillanatnyi frekvencia meghatározására irányult, legközelebb pedig egy másikról fogok írni. Kezdeti törekvésem a jel energiájának a felhasználása felé irányult. Abból a feltételezésből indultam ki, hogy a jel energiája két egymás utáni pontban nagyjából ugyanaz, azaz a pillanatnyi energiaváltozás nem lehet óriási: \displaystyle E_{k+1}-E_k\approx 0, vagyis még jobban általánosítva, az teljes energia változása nulla, azaz deriváltja nulla:  \displaystyle \frac{dE_t}{dt}=0. Felírva a teljes energiát a kinetikus és a potenciális energia összegeként kaptam a következő összefüggést:  \displaystyle \frac{dE_t}{dt}=\frac{2kc\dot{c}}{2}+\frac{2m\ddot{c}\dot{c}}{2}= 0, ahol a  \displaystyle k a jel rugalmassági állandója,  \displaystyle m pedig a rezgőmozgást (körmozgást) végző test tömege. Átalakítások után kapjuk, hogy  \displaystyle k\dot{c}c+m\ddot{c}\dot{c}=0\Rightarrow\frac{k}{m}=-\frac{\ddot{c}}{c}. Azonban tudjuk, hogy a rugalmassági állandó és a tömeg aránya nem más, mint a szögsebesség négyzete, ezért  \displaystyle \omega(t)=\sqrt{-\frac{\ddot{c}(t)}{c(t)}}.

A fent levezetett képlet nem helyes bármely  \displaystyle c módusfüggvényre, pontosabban azokra, amelyeknek az amplitúdója viszonylag gyorsan változik időben. Azonban, amint azt az [2]-ben is bemutatják, a képlet alkalmazható konstans amplitúdójú és változó szögsebességű körmozgás esetén, illetve nagyon lassú változású amplitúdó esetén is. Sajnos az elektrofiziológiai jelek nagy része nem veti alá magát ennek a kitételnek, így a módszer nem alkalmazható a céljaimra.

A kudarc ellenére a befektetett idő és energia mégsem mondható elvesztegetettnek, mert így legalább megtanulhattam, hogy

  1. Nem minden olyan egyszerű, mint amilyennek tűnik;
  2. A jelek modellezése az egyik legnehezebb dolog, amivel valaha is szembesültem;
  3. Sikerült olyan elméleti tudásra szert tegyek a pillanatnyi frekvencia kiszámítását illetően, ami nem sikerült csupán cikkek és könyvfejezetek olvasása közben.

A következő bejegyzésben rátérek arra, hogy hogyan próbáltam meghatározni a módusfüggvény kvadratúráját.

grad_red_right

Research does not consist only of getting the expected results after a tiresome work – as it turned out recently.

My research plan contained the goal of studying the existing methods for computing the instantaneous frequency of a mono-component (or intrinsic mode function – IMF). At the beginning of the first trimester I started to work on this task, and after the study of the Hilbert-transform and the empirical AM-FM decomposition, I started to think about a mathematical model that would help to develop an alternative method for computing the instantaneous frequency. In order to achieve this, I considered the IMFs to be of the form (based on [1]):  \displaystyle c(t)=a(t)\sin(2\pi\omega(t)). Because this is a sinusoidal signal, it was in my right to assume that it has a signal pair, which is cosinusoidal (let this be  \displaystyle c_q(t)), and the two of them are projections of a circular motion with variable amplitude and circular velocity. The cosinusoidal signal is also called the quadrature of the IMF, i.e. there is a phase difference of  \displaystyle \pi/2 between the two of them. Up until here, the idea behind my method is the same as the one behind the Hilbert-transform and the empirical AM-FM decomposition: compute in a way or another the quadrature of the IMF, then compute the instantaneous phase which helps to obtain the instantaneous frequency.

This blog post presents one of the dead-ends I met during my work. In the next post I will present another one. First, I tried to make use of the total energy of the signal (i.e. the IMF) and I assumed that the energy must be equal in two consecutive moments, because it surely cannot vary too much (however this too much was hard to be quantified): \displaystyle E_{k+1}-E_k\approx 0. Taking a further step towards generalising, I can say that the variation of the total energy is zero:  \displaystyle \frac{dE_t}{dt}=0. Writing the total energy as the sum of kinetic and potential energies, I obtained: \displaystyle \frac{dE_t}{dt}=\frac{2kc\dot{c}}{2}+\frac{2m\ddot{c}\dot{c}}{2}= 0, where  \displaystyle k is the spring constant of the IMF and  \displaystyle m is the mass of the body that performs the circular motion. After some transforms I obtained:  \displaystyle k\dot{c}c+m\ddot{c}\dot{c}=0\Rightarrow\frac{k}{m}=-\frac{\ddot{c}}{c}. But it is known that the proportion of the spring constant and the mass gives the square of the angular velocity, thus:  \displaystyle \omega(t)=\sqrt{-\frac{\ddot{c}(t)}{c(t)}}.

The above formula is not correct in any cases. Specifically it is not correct for any \displaystyle c IMF which have a relatively quickly varying amplitude. But, as it is shown in [2], the formula can be applied for constant amplitude IMFs and for those which have slowly varying amplitudes. Unfortunately, electrophysiological signals do not obey this criterion, so this method is not useful at all for my purposes.

Albeit this is a failure, the amount of invested time and energy isn’t a total waste, as at least I learned an important lesson:

  1. Not everything is as simply as it seems at first;
  2. Signal modeling is one of the most difficult tasks I had ever met with;
  3. I managed to gain some theoretical knowledge regarding the computation of the instantaneous frequency that I couldn’t have by only reading articles and book chapters.

In the next post I will describe how I tried to determine the quadrature of the IMF.

Referencia/Reference:

  1. Huang, N. E., Wu, Z., Long, S. R., Arnold, K. C., Chen, X., and Blank, K., 2009, On Instantaneous Frequency, Advances in Adaptive Data Analysis, vol. 1, no. 2, pp. 177–229.
  2. Földvári, R., 1994, Generalized instantaneous amplitude and frequency functions and their application for pitch frequency determination, Journal of Circuits, Systems and Computers, vol. 5, no. 2, pp. 145-165

Nem-stacionárius jelek vizsgálata – Examining non-stationary signals

A legutóbbi bejegyzés óta sikerült egységesíteni a tanulmányozott transzformációkról szerzett tapasztalataimat. A legfontosabb következtetés amire jutottam, az, hogy sem a Fourier-transzformáció (FFT), sem a Short-time Fourier-transzformáció (STFT), sem pedig a wavelet-transzformáció (WT) nem képes hűen tükrözni egy nemstacionárius jel tulajdonságait frekvenciatartományban, illetve idő-frekvencia tartományban. Nyilvánvaló előnyeik mellett a transzformációk a következő hátrányokkal rendelkeznek:

  1. FFT: csupán periodikus jelekre ideális, feltételezi, hogy a frekvencia komponensek végig jelen vannak a vizsgált jelben (egy ilyen jel lenne például \displaystyle x(t)=\sin(2\pi50t)+0.3\sin(2\pi5t), amelyben egyszerre van jelen egy 50 Hz-es és egy 5 Hz-es frekvencia komponens) ezért nem nem lokalizálja a frekvenciákat időben.
  2. STFT: egy időablakkal “pásztázza” végig a jelet, annak reményében, hogy lokalizálni tudja az egyes, időben korlátozott előfordulású frekvenciákat; az ablak mérete befolyásolja úgy a frekvencia-, mint az időtartomány felbontását: ha az ablak mérete nagy, a frekvenciatartomány felbontása nagy lesz, az időtartomány felbontásának róvására, míg a kis ablakméret gyenge frekvenciatartománybeli felbontást és jó időlokalizációt jelent [1].
  3. WT: a transzformáció sikere azon áll vagy bukik, hogy sikerül-e helyesen megválasztani a bázis wavelet-et, amelyet aztan kihúzva és összenyomva, illetve tologatva a jelen megkapjuk a jel wavelet-síkbeli ábrázolását, tehát szükséges a jel tulajdonságainak előzetes ismerete [2].

A fent említett problémákra megoldást jelenthetnek a nemlineáris jelekfeldolgozó transzformációk (például a Wigner-Ville eloszlás [3]) vagy jel-adaptív transzformációk, mint a Hilbert-Huang transzformáció [4]. Ez a transzformáció két lépésből áll:

  1. empirikus módusfelbontás (empirical mode decomposition – EMD): a jelet a nullaátmenetek és lokális szélsőértékpontok alapján “általánosított” harmonikusokra bontja, ún. benső módusfüggvényekre (intrinsic mode function – IMF);
  2. a kapott benső módusfüggvények Hilbert transzformációja segítségével kiszámolja azok pillanatnyi frekvenciáját és pillanatnyi amplitúdóját.

A felbontásból származó módusfüggvények olyan jelek, amelyekben egy adott időpillanatban egy adott frekvenciakomponens van jelen. Ebből kifolyólag összehasonlítható egy olyan forgómozgással, amelynek mind a szögsebessége, mind pedig az amplitúdója időben változó. A módusfüggvények pillanatnyi frekvenciájának és amplitúdójának meghatározása után az eredeti jel így kapható vissza: \displaystyle x(t)=\sum_{i=0}^{N}{a_i(t)e^{j\int{\omega_i(t)dt}}}. Ez felfogható a Fourier transzformáció általános formájaként is; a különbséget az időben változó amplitúdó és frekvencia adja.

Vizsgálataim során észrevettem, hogy a Hilbert-transzformációval meghatározott pillanatnyi frekvencia nem folytonos, időnként ugrások szakítják meg (első ábra). Más módszer után nézve az AM-FM felbontást találtam a legmegfelelőbbnek [5]. Ez a módszer egy IMF-et egy amplitúdó- és frekvenciamodulált jel szorzatának tekinti és megpróbálja szétválasztani a két jelet úgy, hogy az IMF-re egy köbös spline-t illeszt, majd elosztja vele az IMF-et és a műveletet addig ismétli, amig az eredmény egy frekvenciamodulált (FM) jel lesz, amelynek a csúcsai a [-1; 1] intervallumban vannak. Ha az eredeti IMF-et elosztjuk a kapott FM jellel, megkapjuk az amplitúdómodulált részt (AM). Az FM rész és a kvadratúrája segítségével kiszámítható az IMF pillanatnyi frekvenciája.
Azonban ez a módszer sem teljesen problémamentes. Amint a második ábra is mutatja, a spline nem mindig illeszkedik pontosan a jelre, előfordul, hogy egyes szakaszokon kisebb a jelnél, ilyenkor az FM rész nem kerül a
[-1; 1] intervallumba és ezáltal a kapott pillanatnyi frekvencia ugrásokat fog tartalmazni, hasonlóan ahhoz az esethez, amikor a Hilbert transzformáció segítségével számoltuk ki.
Ennek a problémának a kiküszöbölésére a következő algoritmust dolgoztuk ki:

  1. megkeressük azokat a szakaszokat, ahol a spline kisebb, mint a jel
  2. szakaszonként megkeressük a jel és a spline közötti legnagyobb távolságot
  3. a kapott pontot hozzáadjuk a spline-hoz
  4. a fenti lépéseket addig ismételjük, amíg a spline végig nagyobb lesz a jelnél és az FM rész a [-1; 1] intervallumba kerül.

A fenti algoritmust alkalmazva a pillanatnyi frekvencia pontosabban meghatározhatóvá válik, elkerülve a Hilbert transzformáció használatával járó műtermékeket.

grad_red_right

Since my last post I managed to unify the experiences I gained about the different transforms. The most important conclusion I have reached is that neither the Fourier transform (FFT), nor the Short-time Fourier transform (STFT), nor the wavelet transform (WT) is able to reflect intrinsic characteristics of a non-stationary signal in the frequency domain or in the time-frequency domain. Beside their obvious advantages, these transformations have the following drawbacks:

  1. FFT: ideal only for periodic signals as it assumes that the frequency components of the signal are present over the whole time span of the signal (such a signal would be \displaystyle x(t)=\sin(2\pi50t)+0.3\sin(2\pi5t), where both the 50 Hz and the 5 Hz components are available at the same time) and this is the reason why this transfrom can’t localize the frequencies in the time domain.
  2. STFT: “scans” the signal with a time window, hoping to localize frequency components that appear only in some well-defined moments in time. The size of the analyzing window affects the resolution in both the time domain and the frequency domain: if the size of the window is big, the resolution of the frequency domain will increase and the resolution of the time domain will be poor and if the size of the window is small, the time domain will have a greater resolution in the detriment of the frequency domain. This phenomenon is closely related to the uncertainty principle known in quantum mechanics [1].
  3. WT: the transformation if based on the choice of so called base wavelets which will then be stretched, compressed and scaled over the signal in order to transform it to the wavelet domain. Failing to choose an adequate base wavelet, the output of the transform will hardly be useful. This means that a priori knowledge about the signal is needed [2].

The problems presented above could have a solution when using non-linear signal processing transforms (such a transform is the Wigner-Ville distribution [3]) or signal adaptive transforms, such as the Hilbert-Huang transform [4]. The latter one is composed of two steps:

  1. the empirical mode decomposition (EMD): decomposes the signal into generalized harmonic components based on the zero crossing points and the local extrema of the signal, these components being called intrinsic mode functions (IMF);
  2. the intrinsic mode functions undergo a Hilbert transformation in order to compute their instantaneous frequency and instantaneous amplitude.

The IMFs are signals that contain only one frequency component at a given moment. This way, they can be compared to a non-uniform circular motion with variable amplitude. After the computation of the instantaneous frequency and amplitude of the IMFs the original signal may be obtained by using the formula: \displaystyle x(t)=\sum_{i=0}^{N}{a_i(t)e^{j\int{\omega_i(t)dt}}}. It can easily be observed that this is a generalization of the Fourier transform, the only difference being the time-variant amplitude and frequency.

During my investigations, I observed that the Hilbert transform based instantaneous frequency is not continuous in time, but sometimes contains spikes (as seen on middle image of the first figure). While looking for another method, I discovered the empirical AM-FM decomposition [5]. This method decomposes an IMF into an amplitude and a frequency modulated signal by fitting a cubic spline on the IMF and then dividing the IMF with the spline. This division is repeated until the reseult is a purely frequency modulated signal laying in the [-1; 1] interval. Dividing the original IMF with this frequency modulated signal (FM part) one would obtain the amplitude modulated part (AM part). Using the quadrature of the FM part, the instantaneous frequency of the IMF can easily be computed. But even this method fails sometimes as the spline may become less than the signal (as can be seen on the bottom image of the second figure), causing the FM part not to converge to the [-1; 1] interval thus introducing spikes in the final instantaneous frequency values. To solve this issue, we propose the following algorithm:

  1. look for the segments where the spline is less than the IMF signal;
  2. localize the maximal distance between the IMF and the spline for each segment;
  3. the identified maximum is added then to the trajectory of the spline;
  4. repeat the previous steps until the FM part will fit entirely in the [-1; 1] interval.

After the application of this algorithm the instantaneous frequency became much more precise than the one computed with the Hilbert transform, lacking the meaningless spikes.

Tüskés pillanatnyi frekcvenciaértékek (középső grafikon); javított pillanatnyi frekvencia (alsó grafikon)

Tüskés pillanatnyi frekcvenciaértékek (középső grafikon); javított pillanatnyi frekvencia (alsó grafikon)

A spline burkoló kisebb a jelnél (alsó grafikon).

A spline burkoló kisebb a jelnél (alsó grafikon).

Referenciák/References

  1. Gröchenig, K., 2003, Uncertainty Principles for Time-Frequency Representations in Advances in Gabor Analysis, Applied and Numerical Harmonic Analysis, pp. 11-30
  2. Tikkanen, P. E., 1999, Nonlinear wavelet and wavelet packet denoising of electrocardiogram signal, Biological Cybernetics, no. 80, pp. 259-267
  3. Yan, Y. S., Poon, C. C. Y., Zhang, Y. T., 2005, Reduction of motion artifact in pulse oximetry by smoothed pseudo Wigner-Ville destribution, Journal of NeuroEngineering and Rehabilitation, vol. 2, no. 3
  4. Huang, N. E. et al, 1998, The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis, Proceedings of the Royal Society London A no. 454, pp. 903-995
  5. Huang, N. E. et al, 2009, On Instantaneous Frequency, Advances in Adaptive Data Analysis, vol. 1, no. 2, pp. 177-229

Az indulás – The start

Üdv a blogomon!

Ezt az online naplót azért indítottam, hogy időnként beszámolhassak a kutatási eredményeimről és egyéb hírekről. Kutatásom címe “Elektrofiziológiai jelek vizsgálata új jelfeldolgozási módszerekkel” és a TÁMOP 4.2.4.A/2-11-1-2012-0001 azonosító számú Nemzeti Kiválóság Program Határon Túli Kiváló Hallgatói Ösztöndíj segítségével valósul meg. A hét elején értesültem róla, hogy megnyertem a pályázatot, ezért a mai bejegyzés csak egy rövid ismertető a kutatásomról.

A kutatás célja kidolgozni egy jelfelbontási eljárást, amely hozzásegít a pillanatnyi frekvenciák kiszámításához és tükrözi az eredeti jel fizikai tulajdonságait a komponensekben, ezáltal egyértelműsítve a kapott pillanatnyi frekvencia értékeket. Felvetődik a kérdés, hogy miért éppen elektrofiziológiai jeleken szeretnék kísérletezni? Azért, mert az orvostudomány előrehaladtával egyre több diagnosztikai eszköz áll a szakemberek rendelkezésére és ezek mindegyike valamilyen formában mért értekeket, jeleket szolgáltat, legyen szó EKG-ról, EEG-ről, EMG-ről, CT-ről vagy MR-ről. A nagy mennyiségű adathalmazok sokrétű tudást és fokozott odafigyelést feltételeznek egy orvos részéről, ezért megkezdődött a diagnosztizálást segítő berendezések kifejlesztése, segítendő az orvosok munkáját. Ezek az eszközök azonban nem helyettesíteni hivatottak az orvost, céljuk inkább az, hogy kiemeljenek bizonyos jelspecifikus tulajdonságokat, amelyek esetleg elkerülhetik a vizsgáló figyelmét. Ezen eszközök egy másik alkalmazási területe a sürgősségi orvoslás lenne, ahol nem szakorvosok látják el a beteget/sebesültet és telefonos kapcsolatban állnak a központtal, ahol egy szakorvos segíti a mentést. Mindezen eszközök mögött egy-egy megbízható, alacsony tévedési rátával műlödő algoritmus kell álljon.

A kutatás jelen fázisában a különböző transzformációk előnyeit, hátrányait és lehetséges módosítását tanulmányozom, úgy mint a Short-time Fourier-transzformáció és a Hilbert-Huang transzformáció. Az ezekhez tartozó következtetések a következő bejegyzésben kerülnek majd bemutatásra.

grad_red_right

Welcome to my blog!

I have started this blog in order to present the progress and results of my research project entitled “The study of new signal processing methods applied to electrophysiological signals” and it is supported by the TÁMOP 4.2.4. A/2-11-1-2012-0001 „National Excellence Program – Elaborating and operating an inland student and researcher personal support system” program. I found out on the Monday of this week that my research proposal was accepted for the support program so this post will only contain a brief description of my project.

The scope of the research is to develop a signal decomposition method, that helps to compute instantaneous frequencies and reflect physical properties of the original signal in its components, making the frequency vales be unequivocal. You may ask what is the reasoning behind choosing electrophysiological signals as the object of the study? The answer is quite simple: today’s medicine is so advanced that there are many diagnostic tools that provide physicians different measurements – such as the ECG, EEG, EMG, CT or MRI. The tremendous amount of data needs profound knowledge and great attention from the examining doctor so the production of equipment capable of making a diagnosis has started recently. These devices are not meant to substitute a doctor, but to help his work and emphasize certain signal characteristics that are eventually difficult to be observed. Another field of application of these devices is emergency medicine and tele-medicine, where there may be unqualified personnel who would use the recommendations and observations of such a device to communicate them with an expert (doctor) who is at a remote location (typically at the call center). These equipment have a powerful, fault-tolerant algorithm in the background, which also have a low false-positive result rate.

At the moment I’m doing a research with respect to the advantages, disadvantages and possible enhancement methods for different mathematical transforms, such as the Short-time Fourier transform and the Hilbert-Huang transform. In the next post I will make a description of the conclusion drawn from the theoretical study.