# Signal processing by means of WaveletTransform

I updated my question to explain what I want. I have the voltage-time curve from the real industial object. This curve was gotten from the digital oscilloscope:

As you can see it has hight-frequency (kHz) interference. We extracted data from the oscilloscope into the Microsoft Excel file and here it is: http://www.fileconvoy.com/dfl.php?id=g3c9169c7c6eb991d999470396248137bf77d4cabe

Using this data we reconstructed the curve (in MathCAD). As an equivalent to time-axis we have n-axis ("n" - number of the point; every point corresponds to certain moment of time; in all we have 4000 points for 8 ms):

Using Fourier Transofrm we got frequency spectrum (amp-freq curve):

Inasmuch as this signal is nonstationary, the usage of Fourier Transform is not appropriate. So we used Fast Wavelet Transform and the result was quite different (and probably correct). This is an amp-point curve which is similar to the recostructed source curve, but there are displayed amplitudes for 4,7 - 5,5 kHz frequencies (for different wavelet coefficients), which appeared exactly in these moments:

And now by some reason there is a need to process the source data in Wolfram Mathematica, also because it is far and away more powerful than MathCAD. But Mathematica is something new for me in the field of signal processing, and now I have no time to study it. There is almost one week left for us to have this task done. The main purpose of the work is to use different wavelet families and obtain the correct frequency spectrum (amp-freq curve).

We are working with power quality and electromagnetic compatibility. This task is an engeneering task.

So here is the diference between datas. Original data is on the top and the interpolated is below.

-
Are you sure this is an amp-freq curve /the last one/ ? I can clearly see the Волт reading on the y-axis /U/, but am not sure that n is an indication of frequency... – Sektor Jan 15 '14 at 14:18
It's my mistake, I'm sorry. I'll correct it. – Ivan Jan 15 '14 at 14:46
Check now, please. – Ivan Jan 15 '14 at 15:06
I have an idea and am trying to work on it right now. I will probably post an update tomorrow, but I have work + finals. – Sektor Jan 15 '14 at 21:10
@Ivan, there are good answers already; just let me add that from an engineering perspective you have to consider other aspects of your problem: the limitations of your device (inside the "real industrial object"?); is this an ongoing issue or a single instance...? And the socio-engineering aspects: who is going to maintain that code, what sort of environment? Do you need to implement it in hardware or software? You certainly can run Mathematica in a RaspberryPi, which is good news; but it all depends on the context of this problem. – carlosayam Jan 15 '14 at 22:14

Once again I repeat that trying to reproduce something like the curve you are trying to get is going to be inaccurate - There are more suitable representations you can use to get the desired frequency spectrum, but you are the one asking the questions :)

data = Import[
"http://www.fileconvoy.com/gf.php?id=geed872d9b8a38dc6999443310.
4661369818fbd5fa1bf3bc&sts=138977900593152145247060b494909aee48bb2a26b595048647"];


f[list_, pos_] := Module[{x = list}, x[[All, pos]] = Sequence[]; x]

fdata = Flatten@f[data, 1];

cwd = ContinuousWaveletTransform[fdata, MorletWavelet[], {6, 20}, SampleRate -> 500000]

WaveletScalogram[cwd, ColorFunction -> "AvocadoColors", ImageSize -> 700]


ListLinePlot[Total[Abs[#]^2] & /@ Reverse@cwd[All, "Values"], PlotRange -> All,
ImageSize -> 700, BaseStyle -> Thick, PlotStyle -> ColorData[19, "ColorList"]]


WaveletScalogram with the scale axis re-adjusted.

freq = (1000/(# MorletWavelet[]["FourierFactor"])) & /@

ticks = Transpose[{Range[Length[freq]], freq}];

WaveletScalogram[cwd, Frame -> True, FrameTicks -> {{ticks, Automatic}, Automatic},
FrameLabel -> {"Time", "Frequency(Hz)"},
ColorFunction -> "SunsetColors", ImageSize -> 700]


Note that you have to change 500 000 to the value of your SampleRate

NB: Bear in mind that the axes are not scaled ! First observe that the number of octaves multiplied by the voices is equal to the max number on the x-axis and then you can use the relationship to scale the y-axis: $$scaledMagnitude = \frac{2 \times Magnitude}{N}$$ where $N$ is the sample size.

To illustrate the method just described, consider the signal

$$Sin[16 π x] + Sin[4 π x]$$

ListLinePlot[Abs@Fourier@Table[N@Sin[16 π x] + N@Sin[4 π x], {x, 0, 32 π, .001}],
PlotRange -> {{0, 900}, {0, 170}}, ImageSize -> 700]


The frequencies are clearly distinguishable.

Now do the following:

cwt = ContinuousWaveletTransform[Table[N@Sin[16 π x] + N@Sin[4 π x], {x, 0, 32 π, .001}],
MorletWavelet[]]

WaveletScalogram[cwt, ImageSize -> 700, ColorFunction -> "SunsetColors"]


And again

ListLinePlot[Total[(Abs[#]^2)] & /@ Reverse@cwt[All, "Values"], ImageSize -> 700,
PlotRange -> All, BaseStyle -> Thick, PlotStyle -> ColorData[19, "ColorList"]]


There it is ... Do not forget to scale the axes accordingly !

-
As I understand, the y-axis is magnitude, but could you explain do we have on the result's x-axis, please? Also on the last picture you have two peaks with too different magnitude, but on the i.stack.imgur.com/bQarM.png those peaks' magnitudes are almost the same. So what's wrong? I just don't get it. And could you give me a bit more explanations about scaling the axes? – Ivan Jan 17 '14 at 21:35
In addition I want to admit, that my signal is 50 Hz sine + some high-frequency iterference. And the magnitude of the biggest high-frequency impulse can be even compared with the original signal amplitude. But I don't see any magnitude peaks on high frequencies. – Ivan Jan 17 '14 at 21:40
X axis - scales. What's wrong is that as I said trying to approximate the magnitude and frequencies this way is actually not the proper way to do so :D Yes, there's an equation in my answer. According to that equation you have to re-label the y axis, because as you can see the axes are data dependent. That's okay, but the true relationship is that described in my answer. You don't see peaks, because first: few samples; second: you are working with scales on the x axis. Transform the scales to pseudo-frequencies and there you have it. – Sektor Jan 17 '14 at 21:46
I know it is a sine wave. I interpolated the data, re-sampled it, computed the wavelet transform and plotted its scalogram. – Sektor Jan 17 '14 at 21:48
So are those magnitude-scale curves for my signal appropriate? And forgive my stupidity, but I don't understand your equation. We have 3999 samples (for 4000 points), that's right? Peak amplitudes of my sine and of the interference is near 150V both. So for the first curve (DiscreteWaveletTransform) we have 150 = Magnitude * (2/3999). Hence magnitude = near 3*10^5. But we don't have such numbers on the y-axis. – Ivan Jan 17 '14 at 22:39

Since you haven't provided any data I define something like:

data =
Table[ Sin[2 Pi t]
+ 0.86 Sin[97 Pi t] Cos[46 Pi t] Sin[39 Pi t] Cos[19 Pi t] Exp[-102 (1/3 - t)^2],
{t, 0.091, 0.519, 1/4095}];
ListLinePlot[data, PlotStyle -> Thick]


Now let's demonstrate how WaveletScalogram depends on the choice of ContinuousWaveletTransform with various wavelets, we choose a few examples with different ColorFunction, appropriate choice strongly depends on a specific purpose of signal processing and the real data one deals with.

GraphicsGrid[
Table[{ Plot[{Re @ #, Im @ #}& @ WaveletPsi[k[[1]], x], {x, -4, 4},
PlotRange -> All, Evaluated -> True, PlotStyle -> Thick],
WaveletScalogram[ ContinuousWaveletTransform[ data, k[[1]], {7, 12},
SampleRate -> 4095],
ColorFunction -> k[[2]]]},
{k, {{ MorletWavelet[],      "BlueGreenYellow"},
{ MexicanHatWavelet[2], "DeepSeaColors"},
{ PaulWavelet[2],       "SolarColors"}}}],
ImageSize -> 930]


These scalograms should be sufficient to start experimenting with your own data.

-
Thank you for your answer. I have a Microsoft Excel file recieved from digital oscilloscope in field conditions. This file contains 4000 points (voltage-time curve) and I don’t even know should I approximate this data, and if I should, how can I done this task. The signal is nonstationary, so I can’t use Fourier Transform to get the correct frequency spectrum. Thus I should use Wavelet Transform to get amplitudes of every frequency component from my voltage-time curve. But I don’t know how. – Ivan Jan 14 '14 at 19:51
I know exactly that my curve contains kHz-range components. Does this information influence on choosing an appropriate wavelet family (Morlet, Gabor, Mexican Hat, etc.)? – Ivan Jan 14 '14 at 19:52
Can you post a link to the data ? – Sektor Jan 14 '14 at 20:24
fileconvoy.com/… – Ivan Jan 14 '14 at 20:30
@Ivan You should try to update your question explaining all you know about your data, all you have tried to do in Mathematica to process it, what you expect to get, etc. Perhaps I'll be able to update my answer in 2 or 3 days since I'm in a business trip. When you give a clear exposition of the issue we'll take a closer look at your question and expectedly be able to provide a more detailed answer. You should know that signal processing capabilities of Mathematica are relatively new so there are rather few experts in this field. Nonetheless we'll try to help. – Artes Jan 15 '14 at 0:15

An alternative; we will discard the time values.

f[list_, pos_] := Module[{x = list}, x[[All, pos]] = Sequence[]; x]

data = Import[
"http://www.fileconvoy.com/gf.php?id=geed872d9b8a38dc6999443310.
4661369818fbd5fa1bf3bc&sts=138977900593152145247060b494909aee48bb2a26b595048647"];

fdata = Flatten@f[data, 1];

cwd = ContinuousWaveletTransform[fdata, GaborWavelet[4], SampleRate -> 1000000]

freq = (1000000/(#*GaborWavelet[4]["FourierFactor"])) & /@
ticks = Transpose[{Range[Length[freq]], freq}];

WaveletScalogram[cwd, ColorFunction -> "AvocadoColors", ImageSize -> 500,
Frame -> True, FrameTicks -> {{ticks, Automatic}, Automatic},
FrameLabel -> {"Time", "Frequency(Hz)"}]


Now that's a clear WaveletScalogram

data = Import[
"http://www.fileconvoy.com/gf.php?id=geed872d9b8a38dc6999443310.
4661369818fbd5fa1bf3bc&sts=138977900593152145247060b494909aee48bb2a26b595048647"];

dwd = DiscreteWaveletTransform[data]

efrac = dwd["EnergyFraction"]

eth[x_, ind_] := If[(ind /. efrac) < 0.01, x*0., x] /; MemberQ[efrac[[All, 1]], ind]
eth[x_, ___] := x

fwd = WaveletMapIndexed[eth, dwd]

ListContourPlot[Abs@Reverse@Partition[Flatten[fwd[All, "Values"]], 4000],
MaxPlotPoints -> 300, ColorFunction -> "AlpineColors"]


ListContourPlot[Abs@Reverse@Partition[Flatten[dwd[{{_}}, "Values"]], 4000],
MaxPlotPoints -> 300, ColorFunction -> "AlpineColors"]


You can see the difference.

Okay, let's have some fun.

First, the data

data = BinaryReadList["http://www.physionet.org/physiobank/database/ptbdb/patient056/
s0196lre.dat"];

pdat = Take[data, 2000];

cwd = ContinuousWaveletTransform[pdat, GaborWavelet[], SampleRate -> 1000]

WaveletScalogram[cwd, ColorFunction -> "FallColors", ImageSize -> 500]


Oh, that heart ... Does not look good ...

Expanding on Artes' work you can do the following:

data = Table[ Cos[4 Pi t] + 0.3 Sin[55 Pi t] Exp[- 86 (1/3 - t)^2], {t, 0, 1, 1/4095}];

cwd = ContinuousWaveletTransform[data, SampleRate -> 4096]

WaveletScalogram[cwd, ColorFunction -> "RoseColors"]


f = cwd["LinearScalogramFunction"]

Plot3D[f[x, y], {x, 0., 0.999756}, {y, 0.299259, 515.371},
ColorFunction -> "DeepSeaColors", ImageSize -> 500]


Or

ContourPlot[f[x, y], {x, 0, 0.999755859375}, {y, 0.2992592856356853, 515.3711319499473},
ColorFunction -> "SunsetColors", PlotPoints -> 200]


-
How about now ? – Sektor Jan 14 '14 at 21:18
I just don't understand. I need to get frequency spectrum of nonstationary signal, which is performed as list of points using Wavelet Transform. – Ivan Jan 14 '14 at 21:29
Look again at the scalogram. – Sektor Jan 14 '14 at 21:30
Thanks a lot. But could you just explain, how can I get something like intechopen.com/source/html/6880/media/image15.jpeg acoustics.org/press/155th/chapmanfig1.gif using Wavelet Transform – Ivan Jan 14 '14 at 21:43
Whoa, that's a whole different story. Notice how the the x-axis is Freq and the y - Amplitude. Well, the axes of the scalogram are time and scales... You can really easily plot the figure, but you have to construct a new data set consisting of values {freq,amp} – Sektor Jan 14 '14 at 21:52