Calculating Affine Transformation matrix from 4 corner points I'm trying to develop a sudoku solver, using image processing, ocr with neural networks and such. I've got the bare solver working, I can successfully interpret numbers and letters from a picture. All of this is done in GNU Octave. Now I am trying to make it more robust, and I've found this question on stack overflow which greatly helped with most of the image processing. But now I am left with the last part - transforming the single "squares" so that they are uniform and (mostly) square.
What I would like to have: A method that from one set of $(X,Y)$ corner coordinates, e.g.:$(73,86), (117,88), (70,127), (115,128)$ would transform the image with interpolation to the set of coordinates, e.g.: $(0,0), (50,0), (0,50), (50,50)$.

I have looked at the functions supplied with octave and there are these possibilities:

imperspectivewarp
imremap

but both of them require, for me to know the affine matrix.

Is there a thing I've overlapped with these functions, or could someone here point me in a direction to develop such a function?

If you already have the pairs of corresponding two points from each coordinate,
you can use the function fitgeotrans
 to calculate the transformation matrix.

Once you have got the transformation matrix, the transformation can be done with imwarp.

Note that you are going to use Homography transformation rather than affine transformation in this case (because you have specified corresponding 4 points, which may cause trapezoidal distortion).

The code would be like this...

% assume that the input image is in the variable inputImage.
nPixelsX = 50; % assume that there is square image with 50 pixels for each axes after transformation.
nPixelsY = 50;

movingPoints = [73,86; ...
                117,88; ...
                70,127; ...
                115,128]; % set image's corresponding points
fixedPoint = [0,0; ...
              nPixelsY,0; ...
              0,nPixelsX; ...
              nPixelsY,nPixelsX]; % set ideal coordinate's corresponding points
transformationType = 'Projective';  % this should be homology transformation
tform = fitgeotrans(movingPoints,fixedPoints,transformationType); % make transformation matrix

resultReference = imref2d([nPixelsY, nPixelsX]);
desiredImage = imwarp(inputImage, tform, 'OutputView', resultReference); % transform


Time scale and Fourier transform Consider the Fourier transform $F(\omega)$ of the function $f(t)$. The magnitude of $F(\omega)$ depends on $\omega$ and thus also depends on the scale of the $t$-axis. For example, when $f_1(t)$ is a box function which is 1 between $[0,1]$ s and $f_2(t)$ is a box function which is 1 between $[0,1]$ ns = $[0,10^{-9}]$ s, then the amplitude of $F_1(\omega)$ is $10^9$ times larger than the amplitude of $F_2(\omega)$. This is because the time scale is reduced by a factor $10^9$ in the second case.
However, no difference would be obtained when you do a numerical Fourier transform. If you have 100 points of the function $f_1(t)$ equally spaced between $[0,10]$ and also 100 points of the function $f_2(t)$ equally spaced between $[0,10\times10^{-9}]$, then numerically the amplitude of the Fourier transform of both functions is the same (which is not correct). Numerically, you only consider the Fourier transform of the "y-values" $f(t_i)$, so the information about the time scale is lost...
You can solve this problem by multiplying the result with $10^{-9}$ in the second case because the Fourier transform of the rectangular function is proportional with $1/\omega$. However, in general, this proportionality is not present...  How to solve this problem for a general function $F(\omega)$? For example, when the Fourier transform is proportional to $e^\omega$ or proportional to $1/\omega^2$ or ...? Is there a numerical approach to take the scale of the time axis into account?

Most of the classical linear transformations (even filtering) may have three main types of scalings:

natural or none: sum or integral does not have an explicit scaling factor (but the inverse may need one)
in amplitude: because of linearity, at least one reference "unit" signal should have a unit amplitude in the transformed domain, for an easy read of graphs. Typically a unit amplitude sine (which is not $L_1$-integrable) should have a one amplitude at its frequency in the Fourier domain
in energy: to have constant energy in the original or the transformed domain.

Each of these options (sometimes, two of them are equal) depends on the purpose. Fourier  has a specific relation to scale (Where in “Discrete Time Signal Processing” (Oppenheim et al.) can I find the scale-change theorem?):
$$  s(\alpha t) \mapsto  \frac{1}{|\alpha|} S \left(  \frac{f}{|\alpha|}\right) $$
Obtaining scale-, shift- or orientation-invariant features is a long-lasting issue in DSP, and is still a concern in machine intelligence or deep learing.
Along the Fourier line, you can look at the Fourier-Mellin transform or scale represensation (by L. Cohen), and the many feature representations like SIFT, SURF, ORB, BRIEF: Image Matching Using SIFT, SURF, BRIEF and
ORB: Performance Comparison for Distorted


what signed 16 bits little endian mean in pcm audio buffer? I'm trying to read data, I want this data in numbers (int double float something i can operate with in order to make some dsp) from a buffer obtained from a pcm configured as S16_LE, (signed 16 bits little endian) but i don't understand how information is organized inside the buffer. I think 16 bits means that signal values are between -32767 and 32768 but nothing more..

Endianness refers to how multi-byte words are organized, supposing that a byte is an 8-bit word. For big-endian, ordering puts  the most significant byte first and the least significant byte last. This is the converse for  little-endians.

An unsigned hexadecimal number 0x2317 needs  two 8-bytes

big-endian order: 0x23 0x17,
little-endian order :  0x17 0x23.


How to calculate variance of noise on a noisy image (Noise+Signal) I have a long exposure image captured at 70ms. (I1) (As attached)
I also have another image captured at 0.7ms with Gain24. (I2) (As attached)
This image has both noise and signal together.
I would like to compute SNR value by calculating the variance of the signal I1 and variance of noise (I2-I1). Since the gain and exposure are different, the pixel intensity of my I2 image is different from I1, such that (I2-I1) gives not only noise but also a considerable amount of signal.
Is there a way to only compute noise variance from I2 to find SNR. Which also means is there a way to find the true signal in I2

I would assume that I1 and I2 are of the same scene/same lighting conditions/same camera position.

In the ideal world images I1 and I2 would have been related by a constant scale factor (+ random noise) so in this case the way to go would be to determine an optimal scale factor that maps I1 into I2 and then estimating noise variance from the difference. This variance would corresponde to the sum of noise variances of your 2 images, but under the assumption that I1 is significantly less noisy this would give a descent estimate

In the real world there is a thing called camera response curve that maps scene irradiance to pixel intensity values and it's non linear. In fact there can be many more things involved, like gamma correction, custom postprocessing that is done by your particular camera imaging pipeline. Undoing all this mess to get back to linear situation can be an incredible headache. There are software packages and methods for camera response calibration/estimation of other pipeline parameters so in principle it's possible to calibrate the camera prior to taking the images, then undo the camera response curve/other processing making everything linear and after that you can estimate the noise just as in 1). Check camera response curve on google for more details

Applying procedure 1) in the uncalibrated situation gives you a solid overestimate of your noise variance so maybe this would be OK for your purpose.

I had situations where I had to estimate noise level of an image given just a single image. Since typical real world images contain a lot of low frequency regions in practice it's sometimes possible to simply process your image with some smoothing filter (linear or for example median) and then analyse the difference between smoothed image and the original one. For example if at least half of your image correspondes to low frequency regions MAD estimate taken from the difference of your image and smoothed version of your image would give you a usable estimate of noise variance.

Under the assumption of unknown camera response curve you can take a statistical approach to this problem and try to estimate a model with more parameters than just a single scale factor that maps intensities of I1 to I2. Unfortunately I am not ready to give you the exact parametric form of such a model, since I don't see your data.

Hope this helps!


Image Registration by Segmentation Image registration algorithms are usually based on point features such as SIFT (Scale-Invariant Feature Transform).

I saw some references to line features, but I was wondering if it would be possible to match image segments instead of points. For example, given source and transformed image:

I can do edge detection, blurring and Watershed Transform on each:

Regrettably, the segmentation turned out to be too different on each image to match individual segments.

I saw some papers on matching shapes and shape descriptors which are invariant to affine transforms, so this area seems to be promising...

Are there any segmentation methods more robust to affine (or even projective) deformations of the image?

MSER (Maximally stable extremal regions) are regions, not points. And they're invariant to affine transformation. But it's not a segmentation method, strictly speaking

Informally speaking, the idea is to find blobs at various thresholds, then select the blobs that have the least change in shape/area over a range of thresholds. These regions should be stable for a large range of grayscale and geometric transformations.


How to Validate Bilateral Filter Implementation? For this may be a X-Y question, I'm providing the full background story:

I'm trying to test for the correctness for a fixed point implementation of bilateral filter written in C.

My current approach is to compare it's output with a double-precision Matlab version.

This helps me capture issues like insufficient precision, pointer out of bound and uninitialized variables and other (stupid) bugs.

Now the C version matches Matlab version with sufficient precision and with no obvious bugs, but I'm not very confident that Matlab version is correct at all, for it is not fully tested or reviewed by anyone.

Having worked as a tester for a while, I have the feeling that be best way, if not the only way, to test for correctness of the implementation of a signal processing algorithm, is to check for it's properties.

For instance, this CDF 9/7 wavelet example uses "perfect reconstruction" property and "N vanishing moments" property as test cases, and it indeed established a fair amount of confidence that his implementation works.

The only property of bilateral filer is "edge preserving", which isn't a strong proof of correctness.

Does bilateral filter has other properties that can be used as proof of correctness of an implementation?

I'm assuming implementation which is defined by two parameters - Color Range STD and Spatial STD.

Few Properties:

If you set the Color Range STD to be high, the result should be very very similar to Gaussian Blur with the Spatial STD.
If you set the Color Range STD to be very near zero, the output image should be with no change (No smoothing at all).
Good image to test things on is Black White rectangles (Like Checkers).


Stability of system with poles inside unit circle - conflict with differential equation I am trying to understand why a system with a single pole inside the unit circle is stable. For example, take a system with one pole at $z=\frac{1}{2}$. The literature says the system is stable. As a Physics major with not that much intuition for pole-zero plots I first tried to solve the differential equation that this transfer function represents:

$\frac{Y(z)}{X(z)}=H(z)=\frac{z}{z-\frac{1}{2}}$

Which, I re-write to find

$Y(z)\cdot\big(z-\frac{1}{2}\big)=X(z)\cdot z$

Which I transformed to the time domain giving

$\dot{y}-\frac{1}{2}y=\dot{x}$

For the impulse response $\dot{x}=\dot{\delta}(t)=0$, so the solution to my equation looks like

$y_{\delta}(t)=c_1e^{\frac{1}{2}t}$

And this looks quite unstable to me!

Of course, I also tried doing the same by doing a long-tail division of the transfer function, so I got

$H(z)=\sum\limits_{n=0}^{\infty}2^{-n}z^{-n}$

The inverse transform of this gives

$y_{\delta}(t)=2^{-t} = \big(\frac{1}{2}\big)^t$

This, I agree, is a stable system.

I've tried looking all over but I haven't found an example of how to go from the differential equation to a stable system. I have seen that some other questions had answers that included something like "look at the characteristic equation, when $\lambda < 1$ then the system is stable. I suppose I would disagree there - since $e^{\lambda t}$ would only be stable for negative $\lambda$.

What am I missing?

What you are missing is that this is about a discrete-time system, because we're talking about poles and zeros in the complex $z$-plane and about poles inside or outside the unit circle. So there is no differential equation, but there is a difference equation:

$$y[n]=\frac12y[n-1]+x[n]\tag{1}$$

The corresponding impulse response is

$$h[n]=\left(\frac12\right)^nu[n]\tag{2}$$

where $u[n]$ is the unit step sequence.

A continuous-time system with a pole at $s=\frac12$ would indeed be unstable, because for stability, all poles of a continuous-time system must be in the left half plane. In discrete time this corresponds to the requirement of all poles being inside the unit circle $|z|=1$.


Linear Constant Coefficient Differential Equations: Zero-Input and Zero-State responses The solution to a linear constant coefficient differential equation of the form $$\sum_{k = 0}^{N} a_k y^{(k)} (t) = \sum_{k = 0}^{M} b_k x^{(k)} (t)$$ can be written as $y(t) = y_{ZI} (t) + y_{ZS} (t)$ (sum of the zero-input response and the zero-state response) or $y(t) = y_h (t) + y_p (t)$ (the sum of the homogeneous solution and the particular solution). The zero-input response is found by solving the homogeneous equation with initial conditions applied to the solution to the homogeneous equation alone, whereas the zero-state response is found by solving the differential equation with zero initial conditions applied to the solution to the full equation (the sum of the homogeneous solution and the particular solution). However, if we split the solution up to a homogeneous solution and a particular solution (the second decomposition above) we apply initial conditions to the sum of the homogeneous and particular solutions, but not to the homogeneous solution alone first, before adding to it the particular solution. What is the reason for that?

I also wanted to know how this relates to AC steady state analysis in circuit theory. When we do steady state analysis, are we obtaining $y_p (t)$ or $y_{ZS} (t)$? Can you explain why? Is the steady-state response the particular solution or the zero-state response? Is it not one or the other in general?

I would be very grateful if you could explain this in terms of the time-domain description of systems, without resorting to s-plane concepts, poles, etc.

The two approaches should return the same solution. They are just two different ways to get to the same place.

In the ZSR/ZIR method, you are solving two different IVPs - they have the same expression on the left-side of the equality, but different initial conditions and different right-hand sides (one has a $0$, the another one has some function of $t$). In the more mathematical non-engineer like approach (homogeneous/particular), you solve just one IVP with both things: the right-hand side of the equation not set to $0$ and the original initial conditions. Thus one can intuitively see the second approach as some kind of superposition of the two IVPs proposed in the first one.

Just in case you don't believe this yet, let me just write down an example (without doing it in detail as that's not the point of this post) to show that both methods work. Let's solve the following differential equation:

$$\left\{
 \begin{array}{ll}
  y''(t)+4y'(t)+3y(t)=e^{-2t} \\
  y(0)=1 \ , \ y'(0)=1
 \end{array}
\right.$$

The mathematical approach ($y_p+y_h$)

I propose:

$$y_p(t)=-e^{-2t}$$

Check that it satisfies the differential equation and the initial conditions:

$$y_p''(t)+4y_p'(t)+3y_p(t)=-4e^{-2t}+8e^{-2t}-3e^{-2t}=e^{-2t} \ \implies\mathrm{OK}$$
$$y_p(0)=-e^0=1\ \implies\mathrm{OK}$$
$$y_p'(0)=-2e^0=1\ \implies\mathrm{OK}$$

Therefore we conclude that $y_p(t)$ is a valid particular solution. Then we should solve the homogeneous equation. I'm not going to do it step by step, but you can check thay you get something of the form $$y_h(t)=Ae^{-t}+Be^{-3t}$$

If we impose the initial conditions to $y_p(t)+y_h(t)$, the constants turn out to be $A=5/2$ and $B=-1/2$. So the solution can be written as:

$$y(t)=y_p(t)+y_h(t)=\frac52 e^{-t}-e^{-2t}-\frac12 e^{-3t}$$

The engineer approach ($y_{zs}+y_{zi}$)

As stated at the beginning, this approach requires us to solve two different IVPs. The first one is:

$$\left\{
 \begin{array}{ll}
  y_{zi}''(t)+4y_{zi}'(t)+3y_{zi}(t)=0 \\
  y_{zi}(0)=1 \ , \ y_{zi}'(0)=1
 \end{array}
\right.$$

The solution (the generic one, with the constants $A$ and $B$) to this problem is the same as in the previous case, but with different constants (now we only have the homogeneous solution term). In this case, $A=2$ and $B=-1$. Therefore:

$$y_{zi}(t)=2e^{-t}-e^{-3t}$$

Now we have the other IVP:

$$\left\{
 \begin{array}{ll}
  y_{zs}''(t)+4y_{zs}'(t)+3y_{zs}(t)=e^{-2t} \\
  y_{zs}(0)=0 \ , \ y_{zs}'(0)=0
 \end{array}
\right.$$

Note that his is exactly the same IVP as the one stated at the beginning, but with different initial conditions. So we can use the homogeneous/particular approach to solve this IVP as well, and find the constants $A$ and $B$ that suit the initial conditions in this case. We find that

$$y_{zs}(t)=\frac12 e^{-t}-2e^{-2t}+\frac12 e^{-3t}$$

If we sum the two solutions:

$$y(t)=y_{zs}(t)+y_{zi}(t)=\frac52 e^{-t}-e^{-2t}-\frac12 e^{-3t}$$

We got the same result.

Conclusions

As can be seen from the example, both methods arrive to the same solution because they are just two different ways of decomposing the problem into two easier ones, but they are equivalent.

The approach using ZSR and ZIR is more practical in real life. Take for example a circuit. We usually know the initial conditions with no input. When we turn the input on (i.e. we make a sudden change in the circuit), voltages and currents could change abruptly. The ZIR uses these intitial conditions with no input, and the ZSR uses the input when we know all initial conditions are zero (for example, we discharge all capacitors before turning the circuit on). Performing the mathematical approach would not be simple, as we would have to know the initial conditions just when the input appeared, and that can be rather hard to measure or observe. With ZSR/ZIR, we can measure the two things independently and then add them up to get the total solution.

Regarding the steady state question, note that when you solve the differential equation you get the signal $y(t)$. Namely, you can evaluate the output at any instant $t$ you want. Thus, you don't have just the transient response or the steady state one: you have both. You can observe whatever you want as you have the signal as a function of time. The steady state will be reached at $t\to\infty$ (in real life, after some time, when all components in the circuit are acting as the input had been turned on for a long time).


Separating near points I'm programming in matlab. I have processed image:

How can I separate them some thing like this:

In other words, I want to separate points that are near together without morphological operations.

I note that I do not want to use morphological operation such as dilation or erosion.

I have a simple solution for your problem, hope it helps.
You might use a clustering method, as I explain in following:

Assume a certain window size (like $k*k$) as neighborhood.
for each and every pixel count number of bright (non-zero) pixels in the neighborhood.
Choose local maximums (pixels with most populated neighborhood), if there is other local maximums in that region, simply ignore them.
Select that pixel and its neighbors as a cluster.

You might add some modifications to above simple procedure, like merging clusters which are too near to each other and indistinguishable.


DSP Hardware for Detecting an Audio "Click" I am working on isolating a specific audio click in a factory type setting. After my initial analysis I determined that a 2-pole High Pass Filter set between 500 and 800Hz would work best for me.

The timing can be random on when the clicks occur due to human interaction. The shape of the click is a "spike" shape. I do not have the images with me on this workstation but I can get them next week if it helps. Thus the process is not completely automated as a user is inserting the connectors together however the goal is to limit human error so as little interaction with this implementation as possible is preferred.

There will always be 3 clicks. Click A and Click B are very pronounced but Click C is very quiet and up until now a very expensive camera system was used to capture Click C visually.

Potential Microphone: link

However I was curious as to what DSP hardware would be recommended for analyzing such waveforms in a fast and efficient manner.

Side Note: I am leaning towards IIR due to ease of use and tweaking if that is an incorrect assumption please correct me.

Probably the easiest (and most accurate) way to do this is the matched filter approach suggested above.  Basically, you can record each click independently and use those as matched filters against the microphone signal.  If the clicks are always the same, click A will ring up very strong on the A matched filter.  Same with B and C.

To implement the matched filter, you can get a development kit with a processor or DSP chip on it.  It will take some coding to implement the solution (probably C).  In college, I worked on the "TMS320C6713 DSP Starter Kit (DSK)" but you can try Arduino or Beagleboard.  Your system will probably not be highly complex.  As for the microphone, you can get a relatively simple one that connects to the 3.5mm jack on the development kit.  (http://www.amazon.com/Pyle-Pro-PDMIC58-Professional-Handheld-Microphone/dp/B003GEBGA0/ref=sr_1_4?s=electronics&ie=UTF8&qid=1440204270&sr=1-4&keywords=microphone+3.5mm)

A matched filter works by multiplying the incoming signal by the expected stream of data.

Example:
If click A lasts 1mS and your sampling rate is 44.1kHz, you would have a 44 sample matched filter (the click may be much shorter).  You can store this in the code and point multiply your incoming signal by this matched filter.

MATLAB example code:
https://stackoverflow.com/a/19471009/3920284

The maximum value you get when you multiply by your input stream will correspond to the time when your click happened.  You can set a threshold like the example does.

One thing that may change over time is your clicks' sound: as the hardware ages, the clicks will sound different since wearing out may occur.


Why must outputs be the same for a pair of inputs for a causal system? I'm trying to get a good understanding of causality and I can't wrap my head around the first definition of it, as noted in this Wikipedia article:

  Definition 1: A system mapping $x$ to $y$ is causal if and only if,
  for any pair of input signals $x_1(t)$ and $x_2(t)$ such that

  $$x_{1}(t) = x_{2}(t), \quad \forall \ t \le t_{0},$$

  the corresponding outputs satisfy

  $$y_{1}(t) = y_{2}(t), \quad \forall \ t \le t_{0}.$$

Can anyone help explain this?

So $x_1(t)$ and $x_2(t)$ are equal up to a certain time $t_0$, and diverge after that time. What this definition is saying is that, up to time $t_0$, the system will produce the exact same output, without regard for whether the input is actually $x_1(t)$ or $x_2(t)$.

If the system is non-causal, then it may "know" whether the input is actually $x_1(t)$ or $x_2(t)$, because it can "look ahead" in time beyond $t_0$, and it could modify its output accordingly, even before $t_0$.

As an example, consider a lightbulb connected to a battery and to a switch. The switch is initially on, and the lightbulb is also on. At time $t=0$, you decide whether you'll turn off the switch at time $t=60$ (one minute). The lightbulb, being causal, cannot know in advance what will happen, and will stay on at least up to $t=60$. A non-causal lightbulb will know in advance what will happen and may turn itself off at, say, $t=30<60$.

I hope this makes sense -- let me know if you need more clarification.


Bloch equation: adding gradients I have a bloch simulation that I have implemented in Python. The implementation is very simple as:

import numpy as np

GYRO = 267538030.3797  # in radians / s / T
T1 = 0.3
inv_T1 = 1.0/0.3
T2 = 0.05
inv_T2 = 1.0/0.05
M0 = 1.0

# This function is called at various time steps
def bloch(self, M, B):
    # Solve bloch equation in rotating frame.
    # M: Magnetization vector
    # B: Effective magnetic field
    return GYRO * np.cross(M, B) - [M[0] * inv_T2,
                                    M[1] * inv_T2,
                                    inv_T1 * (M[2] - M0)]

Now, during the excitation, the effective field is given by:
So, the pseudocode is:

\begin{align}
B[0] &= \textrm{rf}(t) \cos(2 \pi \omega_{\rm eff}  \ t) \\
B[1] &= \textrm{rf}(t) \sin(2 \pi \omega_{\rm eff}  \ t) \\
B[2] &= 0
\end{align}

Assuming, that the RF excitation is on-resomance i.e. $\omega_{\rm eff} = 0$, I track the evolution of the magnetization of the spin across time and it works and after relaxing for T1 seconds, we have this curve:

So, the transverse magnetization starts to disappear and longitudinal magnetization recovers.

Now, I wanted to simulate a readout gradient which is applied just after the RF pulse. Now, the $X$ and $Y$ values for the effective magnetization is 0 as the RF pulse is turned off. The $Z$ component of the effective field at every time step is given by the product of gradient amplitude and the position. So,

\begin{align}
B[0] &= B[1] = 0 \\
B[2] &= g * p
\end{align}

Now, when I use the same equation to solve this, I get nonsensical results, which is shown by the curve below:

As you can see, the numbers are off the charts and I was expecting an echo at the end. The gradient is a simple trapezoidal gradient and I have checked the units and it all makes sense. I was wondering if my expression for the effective field in the presence of the gradient is correct? I have seen these being represented as rotation matrices as well and I wonder if there is something like that I might have missed.

Update

As per suggestion from M529, I did the following experiments:

# Change the method to ignore relaxation effects
import numpy as np

GYRO = 267538030.3797  # in radians / s / T

# This function is called at various time steps
def bloch(self, M, B):
    # Solve bloch equation in rotating frame.
    # M: Magnetization vector
    # B: Effective magnetic field
    return GYRO * np.cross(M, B)

Now I call the method with the different flip angles (from 0-1200) and note the vector length and it does not maintain the length at larger flip angles. So, perhaps there is a numeric issue. I increased the RF bandwidth for computation reasons. Here is the graph of the magnetization at different flips. It is a weird non-linear effect. It first goes down and then climbs back up uncontrollably.

I am not sure if this is related, but this cross-product is part of the differential equation, something like $\frac{\textrm{d}M}{\textrm{d}t} = \gamma M \times B$... For a true infinitesimal $\textrm{d}t$, this makes sense. For nummerical considerations, this needn't be true. Think of the following:

neglect relaxation for a moment
$\vec{M}_0$ has a certain magnitude (i.e. length)
when you apply an rf pulse, you expect a rotation of $\vec{M}_0$, and a vector's length is preserved under rotations
when you calculate the rotation from the cross product, with several time-steps $\Delta t$, you probably do something like $\vec{M}_{0,n+1} = \vec{M}_{0,n} + \Delta t\;\gamma\vec{M}_{0,n} \times \vec{B}$ (i.e. you compute the change from the previous magnetization and the rf pulse and add it to the previous magnetization vector)
this operation does not preserve the length of the vector: $\vec{M}_{0,n} \times \vec{B}$ is perpendicular to $\vec{M}_{0,n}$. When you add the vectors, the resulting vector is longer than $\vec{M}_{0,n}$ (Pythagorean theorem)

Hence, if you do not take appropriate countermeasures, your magnetization vector will grow. Maybe this is what you see in your simulation: At the end you have an insanely huge signal that does not belong there. You can easily check this via a simple test: Neglect relaxation and apply a huge flip angle, lets say $90000^\circ$ and monitor your magnetization vector length. Certainly it should remain constant. If it doesn't, there might be a bug in your code.

Also try to reformulate your calculations with rotations matrices and check if the results get better.


Symbol rate offset estimation I am confused between symbol rate offset, symbol timing offset and sampling clock offset. Can anyone tell be the difference between these terms.

I'm not sure if there are clear definitions for each term but here are my suggestions:

Symbol rate offset is the difference between the symbol rate entering a receiver and the symbol rate used by the demodulator. It could be measured in Hertz or parts per million (ppm).

Symbol timing offset is the difference between the time that a symbol should have been sampled and when it was sampled. It could be measured in seconds or fractions of the nominal baud duration.

Sampling clock offset is symbol rate offset multiplied by samples per symbol.

John


Acoustical Design considerations for Active Noise Cancellation Project for Air vent I am working on ''Active noise cancellation for air vent'' project. I am
using NI myRIO1900, Audio Technica ATR 3350 omni directional microphones
and Tevion Loudspeakers. I have done the programming using LabVIEW and I
have used FPGA module.

I am using sinusoidal noise source of (200-500 Hz) as test noise signal.

I have implemented LMS on LabVIEW, as a preliminary results I am getting
around 8dB noise reduction but these results are pretty random in a sense.

Only at certain volume levels( if I adjust very carefully), the noise
cancellation is happening and the variation in step size is not giving any
decent change ( again I have to make volume level changes to see some
results).

Over all, I am getting some results but I don't have control over the
results.

This is my set up

[the pipe length is 1 meter]

I know I have to consider secondary path modelling and feedback modelling.
I have actually completed the coding of FBFxLMS also but I don't have any reasonable results with this also.

But my question is how to improve the acoustical design of my setup.

What are the problems using just cardboard material to simulate the environment of air vent?

I have read in 2-3 places (here in 2nd answer) that without proper acoustical design it is not possible to reduce the noise effectively.

Quoting from a paper

  If the acoustical design of the system is not optimized,the digital controller may not be able to attenuate the undesired noise adequately.

The only point I am aware of is using low frequency noises (<= 500Hz) to
consider the sound wave propagation as plane acoustic wave. I am
completely unaware of any other details on how to better my acoustical
design.

I am wondering how to find the limitations of how much a DSP will affect and how much the actual acoustic design. If anyone can provide some reasonable resources to go through, it will be highly appreciated.

I first wasn't sure that this was in the scope of a digital signal processing site, but then decided that it is about limits and environments for doing DSP, so it's quite on-topic.

Your paper (what paper? Please, please, get used to always correctly cite!) hit the spot:

  If the acoustical design of the system is not optimized,the digital controller may not be able to attenuate the undesired noise adequately.

I think you understand quite well what they're saying: You've got to deal with signals in the real world, i.e. with actual, physical, pressure waves. You're doing that by modeling the real world mathematically.

This can of course only be a success if your reality is "easy" enough to fit your model. For example, your ANC probably assumes that there's a position for the noise emitter, and position(s) for the cancellation signal emitter(s). What it probably doesn't involve is that with a bad acoustic design, noise might be leaking out of the back of your emitter, bounce off a distant wall, and come back with a huge, random phase shift.

This is science. In science, you try to reduce the amount of uncontrollable effects as far as possible to make as many things certain as possible.

Your design is a good one, in principle: Cardboard is pretty effective at absorbing sound. However, you weren't really too careful when constructing your system:

gaps. Close them.
Walls. Make them thicker.
backs of speakers, open for absorption and emission of ambient noise. Pad them.
Not too sure about the bottom of your "tube": Is it dampening cardboard or actually the table, which probably is much better at transporting power of a speaker lying flat atop than air is?
Overall system exposed to ambient noise. I'd get a big cardboard, put down a layer of insulation (e.g. styrofoam), put everything inside, fill things up with insulating material.
Is omni-directional really the pattern you're looking for in microphones? This all sounds like you want your noise source mic to be highly directional to only capture the noise, and not the effect your own emissions have at the noise source.

Then, make sure all the components, especially your sensors, work very well:

are your microphones (especially that right in front of the noise emitter) perhaps clipping?
if not 1., then, how clean is the signal your speaker produces? Maybe you want to try and produce a few tones at different volumes, and look at the spectra (coming from an RF background, I'd propose full correlation-based channel sounding, but I don't know how well that works).
Verify your model: unplug your cancellation speakers. Run your device, and record what they would have produced as sound to a file, as well as what your microphones heard. Simulate what the mics at the places where you want your noise to be cancelled should have heard based on these files. Do exactly the same with the  speakers plugged in. Compare simulation with actual microphone recording.


Magnitude spectrum for discrete-time signals I am doing the following question (it is not homework, I am preparing for an exam).

  The time discrete signals $x_1(n)$ and $x_2(n)$ is created by sampling the continuous signal $x_a(t) = \cos (2 \pi 300t) + \cos (2 \pi 600 t)$ with sampling frequency $F_s = 1000\textrm{ Hz}$, $x_1$ with anti aliasing filter and $x_2$ without aliasing filter.

  a) Calculate $X_a(f)$ for the time continuous signal $x_a(t)$.
  b) Calculate the magnitude spectrum for the time discrete signals $x_1(n)$ and $x_2(n)$ i.e $\lvert X_1(f)\rvert$ and $\lvert X_2(f)\rvert$, with $-1<f<1$.

I have done a), I used the table for Fourier Transform to conclude that $$X_a(f) = \frac{1}{2}\Big(\delta(f+300) + \delta(f-300) + \delta(f+600) + \delta(f-600)\Big).$$

But I can not solve b), I believe that I can use the information/answer from a) that is why I provided it.

This looks like a HW problem, so I can only give you a hint. Recall how the sampling process works - it's a convolution with a spike train i.e. a sequence of Dirac delta functions spaced at $F_s$. So the analog spectrum (that you derived in part (a) of the problem) gets replicated many many times centered around integer multiples of $F_s$. Also be careful about aliasing - if the largest frequency in your analog signal exceeds $F_s/2$ it means that this replication process will leak signal energy into the lower frequencies (draw a picture and convince yourself of this).
Finally, analog frequencies in the interval $[-F_s/2, F_s/2]$ get mapped to digital frequencies in the interval $[-1/2, 1/2]$.


Zero Crossing of a Noisy Sine Wave I'm trying to find the zero crossings of a sine wave to transform the sine wave into a square wave. The only problem is that the sine wave is noisy, so I'm getting lots of jitter and false zero crossings.

Can anyone recommend any simple psuedocode or relevant materials? So far I have something like this:

 if (sample[i]>0 && sample[i+1]<0) || (sample[i]<0 && sample[i+1]>0)

Can anyone recommend a more robust method?

You could try low-pass filtering the input signal to get smoother zero-crossings (or even band-pass filtering if you have a good idea of the frequency location of the sine wave). The risk is that if sample-accurate phase information is essential to your application, the additional lag from the filter might be a problem.

Another approach: instead of trying to transform the sine wave into a square wave, what about getting an independent square wave oscillator to align itself in phase/frequency with the sine wave? This can be done with a phase-locked loop.


Why adaptive filter does not work in my application I got a problem when I was trying to denoise a signal. Actually, it is a simple simulation. The signal is the addition of a step signal (The info I wish to get) and a sinusoidal one (the noise I wish to remove). See below
However, I tried different parameters of using the adaptive filter, it simply cannot filter out the sinusoidal noise from the step signal. See figure like this.

Any suggestions will be greatly appreciated!

Below is the matlab code

clear all
close all

%% walking induced noise
t = [1:5000]*1e-2;
f = 0.1;
WalkNoise = 1*sin(2*pi*f.*t);%+1.5*cos(3*pi*f.*t);
WN = WalkNoise + 0.05*randn(size(WalkNoise));

figure
subplot(3,1,1)
plot(t,WN);
title('Noise');

%% signal
h1 = 14; % height of the signal 1
h2 = 18; % height of the signal 2
L = 5000;  % total length of the signal
bp =2500;  % location of break point
x1 = h1*ones(1,bp);
x2 = h2*ones(1,L-bp);
Sig = [x1,x2];
Sig = Sig + 0.1*randn(size(Sig));

subplot(3,1,2)
plot(t,Sig);
title('Signal');

%% walking-induced-noise + signal
NoisySig = Sig + WN;
subplot(3,1,3)
plot(t,NoisySig); hold on
title('Signal + noise');

%% adaptive filtering
figure
plot(t,NoisySig); hold on
title('Signal + noise');

mu = 0.001;            % LMS step size.
ha = adaptfilt.lms(20,mu);
[y,e] = filter(ha,WN,NoisySig);
plot(t,e,'r');
legend('Signal+noise','Filtered using Adaptive filter');

I tried your code change adaptfilt.lms to adaptfilt.nlms
also decrease the step size to 0.0001
These conditions gave me better results.
nlms is better than lms as there is stability in learning filter coefficeints.The lms algorithm could change the filter coefficients drastically.


Kalman filter in practice I have read the description of the Kalman filter, but I am not clear on how it comes together in practice. It appears to be primarily targeted at mechanical or electrical systems since it wants linear state transitions and that it is not useful for anomaly detection or locating state transitions for the same reason (it wants linear state transitions), is that correct? In practice, how does one typically find the components that are expected to be known in advance to use a Kalman filter. I have listed the components, please correct me if my understanding of what needs to be known in advance is incorrect.

I believe these do not need to be known "in advance":

Process noise $\mathbf w$
Observation noise $\mathbf v$
Actual state $\mathbf x$ (this is what the Kalman filter tries to estimate)

I believe these need to be known "in advance" to use a Kalman filter:

The linear state transition model which we apply to $\mathbf x$ (we need to know this in advance, so our states must be governed by known laws, i.e. the Kalman filter is useful for correcting measurements when the transition from one state to another is well understood and deterministic up to a bit of noise - it is not an anomaly finder or a tool to find random state changes)
Control vector $\mathbf u$
Control input model which is applied to control vector $\mathbf u$ (we need to know this in advance, so to use a Kalman filter we also need to know in advance how our controls values affect the model, up to at most some gaussian noise, and the effect needs to be linear)
Covariance $\mathbf Q$ of the process noise (which appears to be time dependent in the wikipedia article, i.e. it depends on the time $k$) - it appears we need to know this in advance and over time, I assume in practice it is taken as being constant?
A (linear) observation model $\mathbf H$
Covariance $\mathbf R$ (which appears to also be time dependent in the wikipedia article) - similar issues to $\mathbf Q$

P.S. And yes I know many of these depend on time, I just dropped all the subscript clutter. Feel free to imagine small letter $k$ to the right and down from each variable name if you would like to.

For some context, let's go back to the Kalman Filter equations:

$\mathbf{x}(k+1) = \mathbf{F}(k) \mathbf{x}(k) + \mathbf{G}(k) \mathbf{u}(k) + \mathbf{w}(k) \\
 \mathbf{z}(k) = \mathbf{H}(k) \mathbf{x}(k) + \mathbf{v}(k)$.

In short, for a plain vanilla KF:

$\mathbf{F}(k)$ must be fully defined. This comes straight from the differential equations of the system. If not, you have a dual estimation problem (i.e. estimate both the state and the system model).  If you don't have differential equations of the system, then a KF isn't for you!

$\mathbf{x}(k)$ is, by definition, unknowable. After all, if you knew it, it wouldn't be an estimation problem!

The control vector $\mathbf{u}(k)$ must be fully defined. Without additional system modelling, the only uncertainty on the control vector may be AWGN, which may be incorporated into the process noise.   Known matrix $\mathbf{G}(k)$ relates the control input to the states - for example, how aileron movement affects the roll of an aircraft. This is mathematically modelled as part of the KF development.

The system process noise $\mathbf{w}(k)$ is also, by definition, unknowable (since it is random noise!). However, the statistics of the noise must be known which, for a plain vanilla KF, must be zero mean AWGN with known covariance $\mathbf{Q}(k)$.  Sometimes, the covariance of the noise may change between samples, but in many cases it is fixed and therefore $\mathbf{Q}$ is a constant.  In some instances, this will be known, but in many instances, this will be "tuned" during system development.

Observations are a similar story. The matrix that relates your measurements to the states $\mathbf{H}(k)$ must be fully defined. Your measurements $\mathbf{z}(k)$ are also known because that's the reading from your sensors!

The sensor measurements, however, are corrupted by AWGN  $\mathbf{v}(k)$, which, being random noise, is by definition unknown. The statistics of the noise must be known, which is zero mean with covariance $\mathbf{R}(k)$. Once again, the covariance may change with time, but for many applications, it is a fixed value.   Often, your sensors will have known noise characteristics from the datasheet.  Otherwise, it's not too hard to determine the mean and variance of your sensors that you need to use.  Yes, this can also be "tuned" empirically.

There are a huge number of "tricks" that can be done to work around the restrictions in a plain vanilla KF, but these are far beyond the scope of this question.

Afterthought:

Whilst googling for "Kalman Filter" results in a million hits, there are a couple of things that I think are worth looking at.  The wikipedia page is a too cluttered to learn from effectively :(

On AVR Freaks, there is an "equation free" intro to the Kalman Filter that I wrote some time ago to try to introduce where it is used for real.

If you're not afraid of maths, there are several books worth reading that are at the senior undergraduate/early postgraduate level. Try either Brown and Hwang which includes all the theory and plenty of example systems. The other that comes highly recommended but I have not read is Gelb, which has the distinct advantage of being cheap!


Damped spring mass system - parameter estimation I've tried to calculate the parameters of a damped spring mass system of the form

$m~ y''(t)+d~y'(t)+c~y(t)=F(t)$

but I have some problems determining the mass m of the system.

The damped spring mass system is given by a MATLAB function

y_out=mass_spring_system(F,Td),

where F is the input force vector and Td the sampling period. The output of the function is already corrupted by noise.

First, I considered the steady state solution for an input force of

$ F(t) = B\cdot \cos(\omega t)$

$\Rightarrow y_p(t) = \frac{B}{\sqrt{(c-m \omega)^2+d^2\omega^2}}\cos(\omega t-\theta)$.

I get the largest response for $\omega=\omega_0=\sqrt{\frac{c}{m}}$, which is the natural resonance frequency.
I estimated this parameter the following using the Power Spectral Density:

NDFT = 2^10;
h = hann(NDFT);
[Pxx,W] = pwelch(F_n,h,NDFT/2,NDFT);
[Pyy,W] = pwelch(y_mean,h,NDFT/2,NDFT);

Ha = sqrt(Pyy./Pxx);
Ha_dB = 20*log10(abs(Ha));

where F_n is white Gaussian noise and y_mean is the ensemble average for 500 realizations. The maximum value should correspond to resonance at the natural frequency $\omega_0$, right?

The second step is the determination of the homogeneous solution:

$y(t) = C\cdot e^{-\gamma t} \cos(\omega_d t + \varphi)$

For the homogenous solution I created an input force vector of

F = [30*ones(1000,1);zeros(1000,1)];

I estimated $\omega_d$ by taking the FFT after averaging the output $y(t)$.
The parameter $\gamma$ I get from determining the envelope of the damped sinusoidal signal using Hilbert transform, taking the log and performing a least square fit using polyfit.

Using the formulas
$\gamma = \frac{d}{2m}$
$\omega_d = \sqrt{\omega_0^2-\gamma^2}$
$\omega_0 = \sqrt{\frac{c}{m}}$

I can calculate the parameters $\frac{d}{m}$ and $\frac{c}{m}$. But how do I get a value for the mass m?

I think I got the answer: I look at the step response of the system.

$H(s) = \frac{1}{ms^2+ds+c}$

The limit of the step response equals my parameter $1/c$

$y(\infty)=\lim_{s\rightarrow 0} s\cdot \frac{1}{s} \cdot H(s) = \frac{1}{c}$

Using $c$ and $c/m$ from the measurements above I can calculate the mass $m$.

I hope this is right, but the results look reasonable...


IDFT of $Y[k]=2X[k]$ for even $k$ If the 16-point DFTs of $x[n]$ and $y[n]$ are given as
$Y[k]=\begin{cases}2X[k], & k=0,2,4,...,14 \\ 0, & k=1,3,5,...,15\end{cases}$, where $x[n],y[n]=0, \forall n<0, n>15$, how can I write $y[n]$ in terms of $x[n]$?

If I take $H[k]=1+(-1)^k$ and multiply it with $X[k]$, I should get $Y[k]$. But when I take its IDFT, it turns out to be zero. How to solve this?

Below is a block-diagram of your signal processing chain (ignoring multiplication by 2).

$$ x[n] \longrightarrow \boxed{ \text{ DFT-N } } \longrightarrow X[k] \longrightarrow \boxed{ \downarrow 2 } \longrightarrow \boxed{ \uparrow 2 } \longrightarrow Y[k] \longrightarrow \boxed{ \text{ IDFT-N } } \longrightarrow y[n]$$

Where sequence $x[n]$ has a length of $N=16$ samples, $X[k]$ is its $N$-point DFT sequence. And $Y[k]$ is related to $X[k]$ as given by
$$ Y[k] = \begin{cases} { 2 X[k] ~~~ , ~~~ k = 0,2,...,N-2 \\ 0 ~~~ , ~~~ \text{ otherwise } \\} \end{cases} $$

I will solve this problem using the duality argument. Eventhough it's not necessary, it can reduce the mathematical burden via bypassing explicit computations.

The duality in DFT stems from the fact that the forward and inverse DFT transforms are (almost) identical, designated with the following:  ($n$ time , $k$ frequency)

$$ x_t[n] \leftrightarrow X_f[k] \iff X_f[n] \leftrightarrow N x_t[-k] $$

We can use the duality property to observe that if an operation  $T$ is performed on the time-domain signals $x[n] \xrightarrow{T} y[n]$ with a corresponding reflection on the frequency domain as $X[k] \xrightarrow{M} Y[k]$, where $M$ is the frequency-domain reflection operator of the time-domain operator $T$ , then the same operation $T$ performed on the freqency domain signals will have the same reflection operator $M$ on the corresponding time domain signals, with a possible scale factor if different sample rates are used on the left and right nodes.

Therefore, to use duality, in the above block diagram, I'll replace the frequency domain relation between $Y[k]$ and $X[k]$ with a time domain relation as:

$$v[n] \longrightarrow \boxed{ \downarrow 2 } \longrightarrow w[n] \longrightarrow \boxed{ \uparrow 2 } \longrightarrow z[n] $$

where in the time domain the operation is compression (downsampling) and expansion (upsampling) whose frequency domain DFT relation is readily available as:

$$ W[k] = 0.5 ( V[k] + V[k+N/2] ) ~~~, k = 0,1,..,N/2-1 $$

where $W[k]$ is $N/2$ point DFT of $N/2$ point compressed sequence $w[n]$, and $V[k]$ is $N$ point DFT of $N$ point sequence $v[n]$. Stated verbally, $W[k]$ is obtained from $V[n]$ by adding its second half onto its first half.

Then
$$Z[k] = W[ (k)_{N/2}] ~~~, k = 0,1,2,...,N-1 $$

Where $Z[k]$ is the $N$ point DFT of length $N$ point sequence $z[n]$. Again stated verbally, $Z[k]$ is merely a repetation of $W[k]$ two times. This is signified by the ramge of $k=0,1,2,...,N-1$ and the modulus operator inside the argument of $W[k]$.

Hence we can find the relation between $Z[k]$ and $V[k]$ as :
$$Z[k] = 0.5 \left( V[ (k)_{N/2}] + V[ (k)_{N/2} + N/2] \right) $$

This relation in the frequency domain was readily available. I didn't use any derivations to achieve it (eventhough at first it was derived by means of DFT manipulations anyway). Now using the duality property of DFT, I assert that, if te time domain signals are replaced by frequency domain signals then their effect on the time domain will similar to the above frequency domain relation, by simply replacing $k$ with $n$. Note that time reversal not needed, as it will be cancelled, but scaling by $2$ is required. So the result will be:

$$ \boxed{ y[n] =   x[ (n)_{N/2}] + x[ (n)_{N/2} + N/2] } $$

where $N=16$, $n=0,1,2,...,15$.

Conclusion. It takes more effort to describe a solution more easily obtained via the duality principle.


Why the following discrete type system,variable accumulator, is time invariant? $y[n]=T{x[n]}=\displaystyle\sum_{k=n-n_{0}}^{n+n_{0}} x[k]$

it is some sort of moving summer which computes $n^{\text{th}}$ output sample by adding all samples lying within length $n_{0}$ around some point $n$ on $k$ -axis (where k is dummy variable) so even if we shift input by amount $n_{0}$ the length of interval of summation remain unchanged which is same as if we're shifting output by $n_{0}$ from which it follows system is Time invariant but how do we prove it mathematically .

my work so far

shifting input by $n_{0}$ , i.e, $T{x[n-n_{0}]}=\displaystyle\sum_{k=n-n_{0}}^{n+n_{0}} x[k-n_{0}]$

now,  $k\to k+n_{0}$ then,

LHS$=\displaystyle\sum_{k=n-2n_{0}}^{n} x[k]\neq y[n-n_{0}]$

here i'm getting the answer Time variant. can anyone help ?

You are overloading the symbol $n$ which doesn't mean the same everywhere in your calculations.  Instead, write down what the sequence of input values is and the corresponding sequence of output values when the input is $x$.  DO NOT use the symbol $n$ anywhere ($n_0$ is OK to use).  So, your answer should look like

\begin{align}
y[-2] &= \cdots\\
y[-1] &= \cdots\\
y[0] &= x[-n_0] + x[-n_0+1] + \cdots + x[0] + x[1] + \cdots x[n_0]\\
y[1] &= x[-n_0+1] + \cdots + x[0] + x[1] + \cdots x[n_0+1]
\end{align}
etc.

Now, shift $x$ by $k$ places (say $k=1$ for starters) and call the result $\hat x$. Repeat the above computation for the input $\hat x$ and call the result $\hat y$.  Is $\hat y$ the same as $y$ except shifted by $k$ places? Repeat for $k+1$ etc until you have managed to convince yourself that regardless of what $k$ you choose, shifting $x$ by $k$ places results in $y$ also shifting by $k$ places.


FIR Filter implementation vs Linear convolution implementation I came across a C code for the FIR filter on one of the websites. It is as follows

 void fir(short * y, const short *x, const short *h, int n_out, int n_coefs)
 {
      int n;
      for (n = 0; n < n_out; n++)
      {
          int k, sum = 0;
          for(k = 0; k < n_coefs; k++)
          {
              sum += h[k] * x[n - n_coefs + 1 + k];
          }
          y[n] = sum;
      }
  }

I used the following input in my main function. Considered n_coeffs as 4 and n_out as 7 (4+4-1). And the coefficients are as follows (after padding)

x ={1,2,3,4,0,0,0};
h={1,2,3,4,0,0,0};

And then, I modified the above FIR code as follows and used it in my function

for (i = 0; i < 7; i++)
{
    y[i] = 0;

    for(j = 0; j < 4; j++)
    {
            y[i] = y[i] + h[j] * x[i - 4 + 1 + j];
    }
}

And the obtained output is {4,11,20,30,20,11,4}

But when I perform convolution operation,

 for(i=0;i<7;i++)
 {
    y[i] = 0;

    for(j=0;j<4;j++)
    {
        y[i] = y[i] + h[j] * x [i-j];
    }
  }

The obtained output is {1,4,10,20,25,24,16}.

Ideally, both convolution and FIR Filter should give the same output, right?

And one more thing I observed is the output will be {4,11,20,30,20,11,4} when the h coefficients are time-reversed before padding i.e. when h={4,3,2,1,0,0,0}.

I am confused a bit, where I am going wrong?

I believe that the code is written in such a way that it gives a correct output only when the coefficients are passed in the reverse fashion. Comparatively, the code will be simpler if we pass the coefficients in reverse fashion.


Good way to detect pulse with known width with background noise? I'm dealing with a light sensor tracking a bursting LED. Unfortunately signal is weak and there are other light sources in environment. I've used optical filters on sensor but there is still some noise.

The LED burst with a fixed duration but the distance and time vary. So I'm interested in finding those 5 pixels wide pulses (the height of pulse may be different). The noise is often a square wave with a much lower frequency.

What algorithm can be used to detect those pulse?

Thanks.

For peak detection a nice method is the following: apply a maximal filter to the data and find the places where the filtered data equals to the original one.

A maximal filter is simply sliding through the data and selecting the maximal element from the sliding window. Formally:

$$g_w[x] = \max\left(f[x-w], f[x-w+1], \dots , f[x+w-1], f[x+w]\right)$$

where $f[x]$ is the original data, $w$ is the half of the window size and $g_w[x]$ is the filtered data. Selecting places where $g_w[x] = f[x]$ are either the local maxima of $f$ or plateaus of local maxima. Obviously we don't want to have plateaus, a single point per peak would be desirable.

The first trick we can do is to apply a minimal filter to the data, select the local minima, then choose places where we detected local maximum (based on the maximal filter) but not a local minimum. This will eliminate long plateaus.

The other thing we can do is to select 'up edges' so places where a non local maximum is followed by a maximum.

To make it clear I generated numeric data from the image you presented and wrote a short code. It uses two constants, one is the window size, the other is the minimum jump we require from a peak.

import numpy as np
from matplotlib import pyplot as plt
from skimage import feature
from scipy.ndimage.filters import maximum_filter

# put your data here
#data = np.array([0, 79, 75, 69, 69, 7, ...])

filter_win_size = 12
peak_intensity_threshold = 10

max_data = maximum_filter(data, filter_win_size)
min_data = -maximum_filter(-data, filter_win_size)

# select places where we detect maximum but not minimum -> we dont want long plateaus
peak_mask = np.logical_and(max_data == data, min_data != data)
# select peaks where we have enough elevation
peak_mask = np.logical_and(peak_mask, max_data - min_data > peak_intensity_threshold)
# a trick to convert True to 1, False to -1
peak_mask = peak_mask * 2 - 1
# select only the up edges to eliminate multiple maximas in a single peak
peak_mask = np.correlate(peak_mask, [-1, 1], mode='same') == 2

max_places = np.where(peak_mask)[0]

fig, ax = plt.subplots()
r = range(data.shape[0])
ax.plot(r, data, 'k')
ax.plot(max_places, data[max_places], 'xr')
ax.grid()
ax.axis((0, 409, 0, 130))
plt.show()


Magnitude of the analytic signal It's well known that magnitude of the analytic signal for narrowband signals gives the envelope. For example we can demodulate AM signal by abs(hilbert(s_AM)) in MATLAB. I'm looking for the exact mathematical conditions which specify when we get the envelope by computing $|s(t)|$. By analytic signal, I mean $$s(t) = x(t) + j\hat{x}(t)$$ Where $\hat{x}(t)$ is the Hilbert transform of $x(t)$, $$\hat{x}(t) = \frac{1}{\pi t}*x(t) = \frac{1}{\pi}\int_{-\infty}^{+\infty}\frac{x(\tau)}{t-\tau}d\tau$$Magnitude of $s(t)$ is $\sqrt{x^2(t) + \hat{x}^2(t)}$ and I don't know how to proceed further.

By definition, the magnitude of the analytic signal $s_a(t)$ can always be considered the envelope of the signal in a purely mathematical sense, simply because the real-valued signal is given by
$$s(t)=\textrm{Re}\big\{s_a(t)\big\}=|s_a(t)|\cos\big(\arg\{s_a(t)\}\big)\tag{1}$$
However, this does not mean that $|s_a(t)|$ has necessarily a physical significance as an actual envelope of $s(t)$. It does have physical significance if the term $\cos(\arg\{s_a(t)\})$ in $(1)$ has the function of a carrier. This carrier doesn't need to be sinusoidal but it must be a narrowband signal, i.e., its center frequency must be large compared to its bandwidth.
If you want a strict mathematical definition of the conditions under which $|s_a(t)|$ is an envelope in a physically significant way, then you must define what you mean by 'envelope' and by 'physical significance'. This is not trivial and probably also not very useful.
If your question is just about the feasibility of envelope detection, then that's a very different question, and one that's much easier to answer. If the modulated signal is given by
$$s(t)=m(t)\cos(\omega_ct)\tag{2}$$
then the magnitude of the corresponding analytic signal equals $m(t)$ if

$m(t)$ has a maximum frequency smaller than the carrier frequency $\omega_c$
$m(t)\ge 0$

These conditions are met for standard AM, but not for DSB-SC.

EDIT: For the signal $(2)$ it's straightforward to show that $|s_a(t)|=m(t)$ if above conditions are met. If $m(t)$ is a low pass signal with a maximum frequency smaller than the carrier frequency $\omega_c$, then the analytic signal is given by
$$s_a(t)=m(t)e^{j\omega_ct}\tag{3}$$
with magnitude
$$|s_a(t)|=|m(t)|\tag{4}$$
Clearly, if $m(t)\ge 0$ then $|s_a(t)|=m(t)$.


The Gradient Operator of a Vectorized Image in Matrix Form I have this optimization problem:

$$ \arg \min_{ X \left( i, j \right) } \sum_{i, j} \left\| X \left( i, j \right) - 255 \right\|_{2}^{2} + \lambda \sum_{i, j} \left\| \nabla X \left( i, j \right) - \nabla Y \left( i, j \right) \right\|_{2}^{2} $$

Where $ X $ is the output image and $ Y $ is the input image.

Let's say the input image is $ y $ and the output image is $ x $ (Transform image to vector) then the problem can be rewritten:

$$ \hat{x} = \arg \min_{x} \frac{1}{2} {\left\| x - 255 \cdot \boldsymbol{1} \right\|}_{2}^{2} + \frac{\lambda}{2} {\left\| {D}_{h} \left( x - y \right) \right\|}_{2}^{2} + \frac{\lambda}{2} {\left\| {D}_{v} \left( x - y \right) \right\|}_{2}^{2} $$

Where $D_h$ is the horizontal Derivative Operator, $D_v$ is the vertical Derivative Operator and $1$ is vector of ones.

Then the solution is given by:

$$ \hat{x} = { \left( I + \lambda {D}_{h}^{T} {D}_{h} + \lambda {D}_{v}^{T} {D}_{v} \right) }^{-1} \left( \lambda {D}_{h}^{T} {D}_{h} y + \lambda {D}_{v}^{T} {D}_{v} y + 255 \cdot \boldsymbol{1} \right) $$

My question is given the input $ y $ how to apply $ {D}_{h} $ and $ {D}_{v} $ to this specific equation.

Thanks for your reply.

It is pretty simple to create those Matrices.
The real issue with them is their size which is enormous for real world images.

For small kernels they are sparse which saves the day.
Indeed for the Derivative Operator, which has only 2 elements, they are highly sparse.

I built them in MATLAB using:

mI = im2double(imread(imageFileName));
mI = mI(11:410, 201:600, 1);

% mI = mI(11:20, 201:210, 1);

numRows     = size(mI, 1);
numCols     = size(mI, 2);
numPixels   = numRows * numCols;

mO = zeros([numRows, numCols, length(vParamLambda)]); %<! Output

vDx = [1, -1]; %<! Matrix is doing Correlation, this for Convolution
vDy = [1; -1];

% Sanity Check
mIxRef = conv2(mI, vDx, 'valid');
mIyRef = conv2(mI, vDy, 'valid');

% mIx = reshape(mDx * mI(:), numRows, numCols - 1);
% mIy = reshape(mDy * mI(:), numRows - 1, numCols);

mDh = sparse(numPixels - numCols, numPixels);
mDv = sparse(numPixels - numRows, numPixels);

% mDx = zeros(numPixels - numCols, numPixels);
% mDy = zeros(numPixels - numRows, numPixels);

tic();
colShift = 0;
for ii = 1:(numPixels - numRows)
    if(mod(ii + colShift, numRows) == 0)
        colShift = colShift + 1;
    end
    mDv(ii, ii + colShift) = -1;
    mDv(ii, ii + colShift + 1) = 1;
end

for ii = 1:(numPixels - numCols)
    mDh(ii, ii) = -1;
    mDh(ii, ii + numCols) = 1;
end
toc();

mIx = reshape(mDh * mI(:), numRows, numCols - 1);
mIy = reshape(mDv * mI(:), numRows - 1, numCols);

mE = mIy - mIyRef;
disp(['Maximum Absolute Error Between Matrix Form to Convolution (Vertical Derivative) - ', num2str(max(abs(mE(:))))]);

mE = mIx - mIxRef;
disp(['Maximum Absolute Error Between Matrix Form to Convolution (Horizontal Derivative) - ', num2str(max(abs(mE(:))))]);

If you load image into mI you see the result is identical to the convolution operator.

The full MATLAB code which implements the solution is at my Signal Processing StackExchange Question 50329 - GitHub Repository.

Remark
In practice usually those kind of equation are solved using Preconditioned Conjugate Gradient Method.
The nice trick is that method requires only the result of the operators applied on a vectorized image. This can be done using the convolution operator itself instead of using the matrix form.


Is there any way to differentiate between impulses? I have a set of signals consisting of impulses at different instant of time i.e. its not a continuous signal. What I want is a mathematically intuitive way to distinguish between any two such signals?

I understand that the time of occurrence of impulses can be used to differentiate between the impulses. But what I need is like a single number that gives a measure of how different two impulse signals are for instance it could be some property of the signal that would be different for any two signals.

Just an additional info the signals are defined on finite duration. I don't know if that info is needed or redundant.

I don't come from a signal processing background so this question could be wrong. Also if I have tagged it wrongly please modify it.

An example of 2 signals would be:

[0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 1 0 0 1 0 1 0 1 1 1 1 1]
[0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 1 0 0 1 0 1 0 1 1 1 1 1 1]

If you are looking at vectors $\mathbf x$ and $\mathbf y$
(presumably of the same length) in which the entries are
mostly $0$s and

with occasional $1$s (unit impulses), a good measure of the difference
between the vectors is the Hamming distance between them, that is, the
number of samples in which $\mathbf x$ and $\mathbf y$ differ
with occasional non-zero entries (impulses of different amplitudes),
then we need to distinguish between impulses of different amplitudes
but the same location. If we use Hamming distance, then the measure
of the difference between
$$\cdots \quad 0\quad 0 \quad 0 \quad 1 \quad 0 \quad 0\quad 0\quad \cdots$$
$$\cdots \quad 0\quad 0 \quad 0 \quad 2 \quad 0 \quad 0\quad 0\quad \cdots$$
and
$$\cdots \quad 0\quad 0 \quad 0 \quad 1 \quad 0 \quad 0\quad 0\quad \cdots$$
$$\cdots \quad 0\quad 0 \quad 0 \quad 9 \quad 0 \quad 0\quad 0\quad \cdots$$
is the same which is not quite capturing our intuition that the two pairs
are more different than the Hamming distance is telling us. Possibilities
to consider in this case are
$$\sum_{i=1}^n |x_i - y_i|,\quad   \sum_{i=1}^n |x_i - y_i|^2,\quad
\mathrm{or}\quad \sqrt{\sum_{i=1}^n |x_i - y_i|^2}.$$
The last quantity is the Euclidean distance between the vectors. Which
measure you want to use depends on how much the measure captures the
features relevant to your application.


Reconstruct under sampled signal I have a signal from a PWM inverter that was sampled at 3,84 kHz. The PWM has a switching frequency of 5 kHz. The pwm signal is feeding an induction motor. If I low-pass filter the PWM signal it should result in a pure sinusoid: A*cos(2*pi*60*t+phi) of 60 Hz, which is the reference signal to the pwm.

But the sampling frequency fs=3.84KHz does not respect the Nyquist theorem that says that the sampling frequency should be twice the highest frequency of the signal, that is, it should be 2*5KHz = 10KHz. So the pwm signal is under-sampled.

By just filtering the sampled signal I believe that aliasing will happen. Right? Then I would not be able to obtain the purely sinusoidal reference.

If I linearly interpolate the PWM sampled signal so that it would go from fs=3.84KHz to fs=10KHz then, I would make the sampling rate equal to the Nyquist frequency. Now, if I filter it, aliasing would not happen anymore. Is this right?

If not, I was wondering if there is any way to filter the pwm signal and obtain the purely sinusoidal reference without aliasing. Interpolation? Resampling? I know a priori that the filtered pwm signal should result in a sinusoid, maybe that could help to reconstruct the sinusoid signal?

Thank you,
L.

PWM is a non-linear process, and it contains rectangular pulses; in theory at least, it cannot be sampled without aliasing because its bandwidth is infinite. In practice, a sampling rate of 5 times the fundamental is enough in many applications; this would mean that you need a sampling rate of at least 25 kHz.

Once the signal is sampled, you can't eliminate aliasing by interpolating. The interpolated signal will still contain all aliases. In fact, it is in general impossible to eliminate aliasing once it occurs.

Even worse, the spectrum of the PWM signal is continuous, which means that the aliases will probably swamp the sideband at 60 Hz. In other words, your problem may have no solution.

EDIT: Would a narrow bandpass filter at 60 Hz work?

Short answer: I don't think it will, because there are likely to be aliases at 60 Hz; these will modify the amplitude and phase of the signal component at 60 Hz.

Longer answer: the spectrum of the undersampled signal is quite complex. First, the analog PWM process is non-linear, which introduces harmonics at multiples of  5 Khz. The actual shape of these harmonics is, I think, not trivial to calculate, but we can assume that they will include sinc-like spectra of very large bandwidth (since the PWM pulses can be very narrow). Adding to this the undersampling at 3.84 kHz, there will probably be complex aliasing happening at 60 Hz.

One thing you can do is look at the analog PWM signal in a spectrum analyzer. Set the tracing to "average" and let it run for a while. Find the bandwdith of each harmonic, and then see if any of them will alias at 60 Hz when undersampled. If (and only if) by a sheer strike of luck none of them do, then the narrow bandpass filter will work as you want.


Some doubts in the matlab code of computing Wigner distribution I found that a discrete form of the Wigner distribution was
$$WD(n,k) = \sum_{m = -N/2}^{N/2}f(n+m)f^*(n-m)e^{-j\frac{2\pi}{N+1}2mk}$$
where $f(n)$ is the signal and time limited within $|n|\le N/2$. I manage to find the matlab code of the Wigner distribution (Sorry that I am not very familiar with the format to put the codes):

 1. function [WD] = EVD_tfrwv(x0)
 2. X=fft(x);
 3. X=[X(1:N/2+1);zeros(N,1);X(N/2+2:N+1)];
 4. x=2*ifft(X);
 5. x=[zeros(N,1);x;zeros(N,1)];
 6. X=zeros(2*N+1);
 7. for k=1:2*N+1   X(:,k)=x(k+(0:2*N));  end
 8. ww=X.*conj(flipud(X));
 9. WW=fft(ww([N+1:2*N+1,1:N],:));
 10. WW=real(WW([N+2:2*N+1,1:N+1],:));

My question is

I understand line 8 tries to get the correlation $f(n+m)f^*(n-m)$ (probably) for different time $n$, but I get confused when line 9 performs the FFT after rearranging the order of these correlation samples. Could anyone tell me why we have to do this?
Line 10 has similar operation.
Why there is a real operator in line 10? I do not see any real operator in the Wigner distribution formula although I understand Wigner distribution is a real distribution.

Thanks for anybody who tries to help!

The FFT comes from the summing operation and multiplication by the
complex pointer term in the formula. Compare it to the formula for
the DFT (or FFT):

$$X_k = \sum_{n=0}^{N-1} x_n\cdot e^{-\frac {2\pi j}{N}kn}$$
This is not similar, real just returns the real part of the array.
The Wigner distribution is real, but because of rounding errors, there is always a very small complex part in the result. This operation removes it and just returns the real part.


Bandwidth of a modulated signal in analog domain If i am using QPSK modulation with data rate 4 Mbps with rolloff factor $0.5$, then the bandwidth will be symbol $\rm rate\times (1 + rolloff\  factor)$. So here it will be $\rm 2\ MHz\times 1.5 = 3\ MHz$
This is baseband data bandwidth.
So should I use analog filter of bandwidth $\rm 3\ MHz$ or should it be more? or should it be more than that?

The analog filter will be more than that such as to minimize further distortion to the passband. The analog filter bandwidth will depend on its purpose (Tx DAC reconstruction, Rx ADC and mixer anti-alias, interference rejection) while the fine passband filtering would be done by the pulse shaping filter in the Tx itself and the matched filter in the Rx, typically both implemented as digital filters.


Sinc Interpolation Artifacts I have written a program that uses sinc interpolation to resample some data. The general algorithm is a that I compute the previous N values and the next N values to get a new sample at a non-integer value.

For example if I have samples:

(0, 8), (1, 5), (2, 3), (3,8), (4,-2), (5, 11)

And I want a sample value at 2.3, I offset the 2.3 incrementally by 1 and take the sinc coefficients and weight them by the sinc function value e.g.:

sinc(-0.7) * 8, sinc(0.3) * 5 + sinc(1.3) * 3 + sinc(2.3) * 8 + sinc(3.3) * -2 + sinc(4.3) * 11

This algorithm fits my understanding of the sinc interpolation concept, and I'm pretty sure there are no bugs.

Unfortunately I'm hear some tin-y artifacts in the sound.

Given the above algorithm is working properly and is a valid way to resample data (a big premise), is there any reason sinc interpolation would create these artifacts? When I use the same method to down sample by a large enough factor the artifacts to go away.

The width of the Sinc (both the main lobe and the window) needs to correspond to a low pass filter with a cutoff frequency suitable for the bandlimit of the data and the target sample rate.  Depending on the number of coefficients, a non-rectangular window may be needed.


Possibility of using Power allocation with Maximum likelihood estimation I have a basic question about maximum likelihood (ML) estimator and its implementation with power allocation.
Suppose we are using $MIMO$ system with $4$x$4$ transmit and receive antennas, so the received signal which includes the signal + noise + ISI can be written :
$y = HX + n$
$H,X, n $ are the channel, transmitted signal and noise, respectively.
suppose that only first antennas is active and other antennas are not active, so the transmitted signal will be as below:

My question is about the possibility to use power allocation to detect the active antenna, which means the estimated signal will be as:

I am asking if we can detect the active antenna in that way, is that possible ? Is there another method we can use instead of that ? Is there any document explain that ?
Thank you.

What you mention could be right in low modulation order, such as QAM but in 16-QAM and 32-QAM, the amplitude is changed also, I mean you will have different amplitude levels, so you can't implement that detector.


Zero-frequency response of an IIR transfer function I have a question about the zero-frequency (DC) response of a IIR filter. The coefficients of the 3rd order IIR filter are as follows,

% polynomial coefficients for the numerator of the IIR transfer function
b = [0.0066324, -0.0130292, 0.0063988];

% polynomial coefficients for the denominator of the IIR transfer function
a = [1, -1.4158855, 0.4158913];

Note that the roots of a are all less than unity, but one is quite close to the unity. I am trying to understand why MATLAB's filter() function produces the following results:

rng default;
x = rand(15,1);
y = filter(b,a,x);
mean(y)/mean(x)

>> 0.0012

In the last statement, I computed the ratio of the input and output mean, MATLAB returns 0.0012.

However, from the zero-frequency of the transfer function $H(z)=0.3448$ (i.e., sum(b)/sum(a)=0.3448), I have expected that mean(y)/mean(x)=0.3448. Checking this,

[H,W] = freqz(b,a);
H(1)

>> 0.3448

Have I misunderstood something here about MATLAB's filter() function, or other misconception of the zero-frequency response of a transfer function? Thanks to anyone who is able to provide an insight to this problem.

Because of the pole near the unit circle the impulse response doesn't really settle until after tens of thousands of samples. In comparison, your input of 15 samples is very short. Try making it longer:

x = rand(15E6,1);
y = filter(b,a,x);
mean(y)/mean(x)

>> 0.34252

Also, the filter starts in a state that corresponds to an all-zeros input history. If you precondition it with a run of the actual input, the result will converge faster:

x = rand(15E6,1);
y = filter(b,a,x);
mean(y(5E6:length(y)))/mean(x(5E6:length(x)))

>> 0.34483


What's the best method to extract vocals from music? I want to achieve something like this https://soundcloud.com/user7421656/sonicworx-isolate-preview .

I've researched it and a few papers mention the use of "descriptors" on the audio and use thresholds on the descriptors to determine singing sources. But before I dive in and go hard out on how these people achieved it, I wanted to make sure it's the best way to achieve an effect similar to the above link.

What's an effective method to achieving singing source separation?

The example sounds like pitch tracking combined with comb filtering. Using method like this could be implemented by "drawing" as well, as in SPEAR analysis and re-synthesis program.

If the vocal is in the middle and other instruments are not, then you can also try center channel isolation:

http://www.virtualdub.org/blog/pivot/entry.php?id=102

The methods however, sound horrible in addition to being inaccurate and full of artifacts. Better solution to the problem would be a program that could effectively guess which part of the signal belongs to the voice and which doesn't. Much more difficult for sure though.


How does the 'phase' of a cosine wave become 'pi' when it transitions from positive to negative? I was watching this video lecture on filters where this professor says that:

  The phase of a cosine wave $\cos(\omega)$ is 0 till $\omega=\pi/2$. As soon as $\omega$ goes beyond $\pi/2$, i.e when $\cos(\omega)$ becomes negative, its phase becomes $\pi$.

I couldn't understand this. How can $\cos(\omega)$ have phase in itself if it is just a plain normal cosine function? What is the physical interpretation of this phase?

The statement in the video is, at best, misleading. Normally, when people speak of the "phase" of a sinusoidal function $\cos \omega$, they are talking about the value of $\omega$. For instance, when we say that a wave undergoes a "phase shift", we mean that the value of $\omega$ has changed. More information on this definition of the word "phase" is at https://en.wikipedia.org/wiki/Phase_(waves).

As noted by Jazzmaniac in comments, we could also try to interpret the statement by regarding the value of $\cos \omega$ as a complex number $z$. In polar representation we can write $z = |z| e^{i \phi}$ where $\phi$ is the phase of the number $z$. This use of the word "phase" is misleading in this context. For instance, applying a "phase shift" would mean that we multiply $\cos \omega$ by some complex number $e^{i \phi'}$. That certainly wouldn't implement a phase shift of the sinusoidal wave, though, in the usual sense.


What is the best input for de-noising autoencoder for sound data? I am currently trying to build an autoencoder to de-noise audio data.
However I have not found any good articles explaining about the input to the autoencoder, i.e. feature vector.
As in speech recognition/text up projects, MFCC is known to perform good as a input data/feature vector to neural network.
I would like to know what kind of data is need for denoise autoencoder for audio/sound data.

Thank you so much for your kind help in advance.

As you mentioned MFCC features are one of the best features to represent audio as it captures both the time and frequency variations in the audio clip.You can get more details about MFCCS features in the below link:
http://practicalcryptography.com/miscellaneous/machine-learning/guide-mel-frequency-cepstral-coefficients-mfccs/

You can import python_speech_features in python and use them to directly find the MFCCs of audio which will be in a matrix form.

You can feed the noisy audio features as input and the clean audio features as output to Deep Denoising AutoEncoder(DDAE) and then you can test on test datasets.

If required you can reconstruct the predicted mfcc output of the DDAE and check the extent to which your architecture works.

You can use Griffin Algorithm for this. Refer the link below for details about the reconstruction.
https://timsainb.github.io/spectrograms-mfccs-and-inversion-in-python.html


Convolution of shifted signal If $y(t) = x(t)*h(t)$, then what is the expression for $y(t+a)$?

Is it $x(t+a)*h(t+a)$ or $x(t+a)*h(t)$?

From your confusion of $x(t+a) \star h(t+a)$ vs $h(t) \star x(t+a)$
I guess that a little help on the argument manipulations on functions and convolutions could be appropriate here, moving in simple examples:

First let us express the usual simplistic case.  Consider the relation:
$$ y(t) = h(t) \cdot x(t) + g(t) \tag{1}$$

then manipulations on the argument $t$ is applied to all functions on both sides

$$ y(t+a) = h(t+a) \cdot x(t+a) + g(t+a)$$

or an arbitrary transform on $t$ would similary be:
$$ y(\phi(t)) = h(\phi(t)) \cdot x(\phi(t)) + g(\phi(t))$$

Now consider the case where two functions convolved to produce the third:
$$y(t) = \int_{-\infty}^{\infty} h(\tau) x(t-\tau) d\tau $$

which is abbreviated as
$$ y(t) = h(t) \star x(t) \tag{2} $$

Now be careful to interpret the case-2. The variable $t$ shows in all functions as an argument but you may not apply the transform on $t$ as you did with the case-1, so assume you have an arbitrary transform on $t$ as $\phi(t)$ then

$$ y(\phi(t)) \neq h(\phi(t)) \star x(\phi(t)) \tag{3} $$

For example, as in your case, if $\phi(t) = t+a$ then you get

$$ y(t+a) \neq h(t+a) \star x(t+a)$$
but
$$ y(t+a) = h(t) \star x(t+a) = h(t+a) \star x(t) $$

The justification of this can (only) be seen when you consider the integral definition of the convolution operator:

$$
\begin{align}
y(t) &= h(t) \star x(t) \\
& = \int_{-\infty}^{\infty} h(\tau) x(t-\tau)d\tau \\
y(t+a) & = \int_{-\infty}^{\infty} h(\tau) x((t+a)-\tau)d\tau \\
& = h(t) \star x(t+a) \\
\end{align}
$$

Note that since the live variable inside the integral only happens in one function ($x(t-\tau)$ in this case) then a change in $t$ will only affect one of them and you get:

$$ y(t+a) = h(t) \star x(t+a) $$
or from commutativity of convolution you get
$$ y(t+a) = h(t+a) \star x(t) $$

So this provides the answer you were looking for. However, it's not over. Because the following case represents an exception:

$$ y(-t) \neq h(t) \star x(-t) $$
but
$$ y(-t) = h(-t) \star x(-t) \tag{4} $$

So how to see this case-4. Again, using the integral definition :

Assuming that $y(t) = h(t) \star x(t)$, then compute the convolution between two new signals $g(t)=h(-t)$ and $z(t)=x(-t)$ as:
$$
\begin{align}
w(t) &= g(t) \star z(t) \\
& = \int_{-\infty}^{\infty} g(\tau) z(t-\tau)d\tau &g(\tau)=h(-\tau),z(t-\tau)=x(-(t-\tau)) \\
& = \int_{-\infty}^{\infty} h(-\tau) x(-(t-\tau))d\tau  &\text{ let } \tau'=-\tau \\
& = -\int_{\infty}^{-\infty} h(\tau') x(-(t+\tau'))d\tau'  &\text{ replace } \tau'  \text{ with } \tau \\
& = \int_{-\infty}^{\infty} h(\tau) x(-t-\tau) d\tau  \\
& = y(-t) \\
\end{align}
$$

hence we conclude that $h(-t) \star x(-t) = w(t) = y(-t) $. As stated before, you must always consult to the (explicit) integral definition to decide on the correct functions used in the convolution operator.


fourier transform of smart-phone accelerometer in matlab I'm new to matlab. I want to process my smart-phone accelerometer data in matlab. I know Matlab let's you connect your phone via USB cable to see accelerometer data in realtime. But according to some reasons I don't want to use this method. Instead I've developed an application that saves accelerometer data. I want this application to produce a ".m" file and use this file to represent x-y-z signals. The problem is that the accelerometer data rate is not constant. for example: "data-1" corresponds to time "0.1", "data-2" corresponds to time "0.2", and "data-3" corresponds to time "0.8". You can see the time interval can change. I need to get fourier transform of these signals, so I can not use the simple matlab vector to represent them. How can I represent these signals? How can I get the correct fourier transform of these signals?

Two methods occur to me.

Perform a non-uniform DFT. You can use Google to find MATLAB implementations of it, or try to write one on your own. This file seems to be legit, although I haven't tried it myself.
If you have the Signal Processing Toolbox, then you can use the resample() function to interpolate your non-uniformly sampled vector linearly onto a vector of uniformly spaced instants, and then perform the usual DFT on the latter.


Signal's Fundamental Frequency is different from Plotted Signal I've been attempting to plot the following function using MATLAB:
$$
x(k)=\sum_{n=11}^{50} \sqrt{n} \sin (2n\pi k) +\sum_{n=1}^{40}\sqrt[3]{n} \sin (3n\pi k)
$$

Note that $k$ is a continuous variable.

The problem I am facing with MATLAB is when I plot it, I get a signal which has a fundamental period of 2 seconds, and not 1.5 seconds.

Have a look below:

I've tried a higher sampling rate:

So, what is going on wrong? Why isn't the fundamental period 1.5Hz?

With $k$ a continuous variable, the frequencies of the sinusoids in the first sum are $f_n=n$, i.e. they are integers. The frequencies in the second sum are integer multiples of $\frac32$. The highest fundamental frequency of which all those frequencies are integer multiples is $f_0=\frac12$. The frequencies in the first sum are given by $2mf_0$ (with integer $m$), and the frequencies in the second sum are $3mf_0$. With a fundamental frequency of $f_0=\frac12$, the fundamental period is $T=2$, as observed by you.


How to properly convert from dB back to amplitude I am new to signal processing.

I have a signal that I want to convert to dB, process and then transform back.

My understanding of the transformation to dB scale is

$X_\mathrm{db} = 20 \cdot \mathrm{log}_{10}(|X|)$ and therefore $|X| = 10^{X_\mathrm{db} / 20}$

However, when calculating the reverse of my audio the sample sounds distorted. See reproducible example below

Can someone explain why this occurs?

import wget
import os
import librosa
import IPython.display

filepath = str(os.getcwd()+"/metallic-drums.wav")
wget.download("https://tinyurl.com/wx9amev", filepath)

MD, sr = librosa.load('metallic-drums.wav')

#Playing the orginal file
IPython.display.Audio(MD, rate = sr)

#Chaning to db
MD_db = librosa.amplitude_to_db(MD)
#Changing back to amp
MD_r = librosa.db_to_amplitude(MD_db)

#Play the result
IPython.display.Audio(MD_r, rate = sr)

You have the abs(X) in the formula, which is not a reversible operation.

All the negative values in the original audio sample has been converted to positive.


Whether Fourier transform formula be considered as Convolution or Correlation? The expression for Fourier transform is given as
$$ F(\omega) = \int\limits_{-\infty}^{+\infty} f(t) \ e^{-j \omega t} \ dt  \tag{1}$$
Now, let one function be $f(x)$ and other be $e^{j\omega t}$ then their convolution at $p(0)$ is given as
$$ p(0) = \int\limits_{-\infty}^{+\infty} f(t) \ e^{-j \omega t} \ dt  \tag{2}$$ same as in $(1)$.
Also, if we split the exponential term into sinusoidal,the expression $(1)$ becomes
$$ F(\omega) = \int\limits_{-\infty}^{+\infty} (f(t) \cos{\omega t} -j f(t)\sin{\omega t})\ dt  \tag{3}$$
Let $f(t)$ and $\cos{\omega t}$ be 2 functions. So their correlation  at value
$q(0)$  is given as
$$ q(0) = \int\limits_{-\infty}^{+\infty} (f(t) \cos{\omega t}) \ dt  \tag{4}$$ same as in $(3)$. (we can do with f(t) and sine function also.)
So I want to know

Whether Fourier transform formula be considered as Convolution or Correlation?

Also, what is meaning of $p(0)$ and $q(0)$ terms ?

Correlation and convolution are basically the same operations. You can express the cross-correlation of two functions $f(t)$ and $g(t)$ by a convolution:

$$R_{fg}(\tau)=f(\tau)\star g^*(-\tau)$$

where $\star$ denotes convolution, and $*$ denotes complex conjugation.

If you evaluate the cross-correlation at $\tau=0$ you get the inner product of $f(t)$ and $g(t)$, and that's exactly what the Fourier transform is: it is the projection of $f(t)$ on the complex exponential $e^{j\omega t}$.

So, to answer your questions:

It is actually an inner product. Due to the equivalence of convolution and correlation, it can be seen as both, evaluated at $\tau=0$.
As already mentioned, the meaning of these terms is an inner product.


How to transform lowpass fir filter to bandpass fir filter without using a built in function in matlab The FIR low-pass filter was designed in MATLAB which characteristics are listed below. Coefficient of this filter was written in variable h. Basis on this filter design a band-pass filter with central frequency 1/5(normalized to fs) keeping the same gain and bandwidth. Give the listing in MATLAB(no using buitl-in function) which allow to set down coefficients of designing filter hx.

Thanks in advance

Mathematically, applying a FIR with impulse response $h_\mathrm{lpf}[n]$ to digital signal is convolution:

$y = x * h_\mathrm{lpf}$, or thanks to the properties of the (discrete) Fourier transform,

$Y = X\cdot H_\mathrm{lpf}$, as convolution becomes multiplication.

Now, making a bandpass out of a low pass can be modeled by shifting the Frequency response $H_\mathrm{lpf}$ in frequency domain. "Shifting" can be represented by a convolution of the low pass filter with a dirac impulse at the desired center frequency:

$H_\mathrm{bpf}= H_\mathrm{lpf} * \delta_{f_\mathrm{center}} $

Again, convolution becomes multiplication when transformed to time domain. The inverse (discrete) Fourier transform of a dirac at ${f_\mathrm{center}}$ is a complex oscillation $e^{2\pi{f_\mathrm{center}}n}$, so this becomes

$y[n] = x * (h_\mathrm{lpf} \cdot e^{2\pi{f_\mathrm{center}}n})$


Continuous Wavelet Transform time vector in python I have a signal sampled at 128 Hz. I used to extract features with the spectrogram function and I decided to upgrade my algorithm and I'm trying to analyze it using  Continuous Wavelet Transform  (pywt.cwt) in python.
this function has only 2 outputs: coefficient and frequency, while spectrogram returns the time vector as well.
Similar to what is said in theory, the resulting frequency vector has varying differences between each two samples (high density at low frequencies, and vice versa). According to the uncertainty principle, the time vector also has to change its density inversely to the frequency vector, is not it? My understanding is based on the famous graph in the picture attached.

If I'm not mistaken, the time vector (that is not returned) is just the original time vector of the signal, meaning the difference between each two adjacent samples is constant.
Please help me understand this mismatch.

Similar question here Which time-frequency coefficients does the Wavelet transform compute?

The picture you've shown is used for DWT such as pywt.wavedec, not CWT.

CWT is a continuous function; it exists at all points in the time-scale plane.

pywt.cwt produces a 2D array with constant numbers of columns and rows, so it's a regular sampling of the CWT, at whatever scales you specify, for every sample of the input:

shape(y)
Out[5]: (512,)

shape(freqs)
Out[6]: (128,)

shape(coef)
Out[7]: (128, 512)

I think the DWT's coefficients contain just enough information to reproduce the original exactly with an inverse DWT, while the CWT output has a lot of redundant information and is maybe missing some information.


Measure object planar displacement on slight moving or vibrating plane I'm working in an application where I'm tracking an object 2d displacement on a given plane and the camera is already calibrated (intrinsic) and the lens distortion has been corrected.

In order to get the object displacement measurements I'm thinking to calculate/find the homography between the two planes (world and image).

I'm assuming that as long the object movement remains within the calculated plane everything will be fine, but now turns out that the plane might move/vibrate slightly (~ +-5mm), would I still getting good results? or is there another better way to do this (calibration) on a moving/vibrating plane?

A quick setup sketch:

...turns out that the plane might move/vibrate slightly (~ +-5mm), would I still getting good results?

I do not see why not. And this is making reasonable assumptions about the equipment.

You have a monocular camera, tracking motion from a viewpoint that is vertical to the direction of motion.

Now, if the plane moves purely in the Z direction (that is, on an axis that is vertical to the lens), then the only thing that will change is the size of the object.

The problem this might create is in the Signal To Noise ratio of the tracking. In other words, if you are tracking a marker against a background and the marker's size "pulsates", then the tracking estimate of the frames where the target is relatively larger will be better (within reason).

By how much will your marker / object size change, depends on the focal length of your lens.

If the plane does not move purely in the Z direction but sort of oscillates back and forth as if it was a pendulum, then it introduces an additional perspective error to the measurement that is proportional to $\sin(\theta)$ where $\theta$ is the angle of oscillation.

In other words, a monocular camera imaging a vertical line segment cannot know if the imaged size is true or the projection of a much larger vertical line segment that stands at some angle $\theta$ towards the lens.

In general however, the "magnification error" (first case) grows smaller with the length of the lens (lower in long lenses) and vice versa. There of course, you have to balance the field of view too so that it includes the variation you are trying to measure.

If the equipment you are using has specific characteristics regarding focus and image formation, you are going to have to take those into account of course.

Hope this helps.


Designing an efficient curve-matching algorithm I'm currently designing a curve-matching algorithm and as I already explored many ideas I'm requesting your help. So, if you got some advices about how to handle this problem, feel free to answer!

The Context

I got a flash LIDAR spatio-temporal response to one short impulsion. The data acquired is stored in a 3-D matrix:

First and second dimensions are spatial ones. As in a 2-D image they represent azimuth and elevation of a pixel.
Third dimension is the temporal one. For each pixel the illumination is recorded for a fixed duration and stored in this dimension. Typically, the curve observed is a Gaussian and its parameters (variance, max amplitude and delay) are affected by the LIDAR surrounding environment.

To be more specific the laser pulse is distorted by the environment. The multi-path of light impacts the Gaussian variance, the object distance to the RADAR impacts the delay and the amplitude and the reflectivity of the object impacts the amplitude of the Gaussian.

The Aim

As in the real world objects are contiguous and homogeneously distributed in space, neighboring pixels are more likely to have a similar reflectivity and distance to the LIDAR thus neighboring Gaussians are very similar.

What I'm trying to do is to quantify the similarity of a pixel impulse response to its neighborhood (3x3 neighborhood for instance). This is very close to the concept of gradient, I want to detect pixels that are parts of an homogeneous objects and ones at the border of an object where there is a disruption.

For now I've implemented two different algorithms. One that is similar to a an image derivation adapted to the 3-D structure, with a 2-D Laplacian Filter. The second algorithm is based on the measure of similarity using Sum of Absolute Differences, also adapted to the data 3-D structure.

The Problem

These algorithms work somewhat fine but since the signal can be noisy and sometimes very far from a Gaussian shape I think and I hope better algorithms can be implemented.

My research tracks are focused on curve-fitting before comparison (Levenberg–Marquardt algorithm for instance) for now, but this proved very computationally time-consuming!

If you know about other methods or concepts that can help in resolving my problem, feel free to participate :)

Thank you!

EDIT:

Matrix size is 32x128x128

The two most obvious things you can try are:

Fitting a Gaussian to your data and then clustering their parameters
Estimate the similarity of waveforms directly and then try to cluster that

Since you know that the return waveform conforms to a Gaussian, it is better to use a method that takes this into account.

So, basically, for every pixel time course, you will end up with a fitted $\mu, \sigma^2$. You can then use these two, plus any other target dependent feature to characterise the surface.

  These algorithms work somewhat fine but since the signal can be noisy and sometimes very far from a Gaussian shape I think and I hope better algorithms can be implemented.

The fact that you know that your data is supposed to conform to a specific model is a great deal of knowledge and you should try to take that into account as much as possible. "...sometimes very far from Gaussian..." doesn't always mean "not Gaussian". You might for example be getting two or more returns in the same time-course which means that you might have to apply something like a kernel density estimator or an iterative fitting process to pick up each individual return. For more information, please see this paper.

The second option is to estimate "similarity" of waveforms directly. In this case, given two time series $U,V \in \mathbb{R}^\mathbb{N}$, you are looking at a function $F(U,V,\Theta) \rightarrow \mathbb{R}$ which simply returns a number that is proportional to how "similar" these two waveforms are, given some additional parameters in $\Theta$. Most of the times, "similar" is another word for "co-varying".

The most obvious of the possible $F$s is the discrete cross-correlation. In that case $f_{xcorr}(U,V)$ returns a number between $-1..1$ which tells you how similar $U,V$ are. Another possible $f$ is Coherence. In this case, parameters include the frequency component of which you want the coherence calculated at (e.g. $f_{coher}(U,V,k)$ (for the $k^{th}$ component) which returns a number between $0..1$. Another possible $f$ is Mutual Information with is a slightly different take on "similarity" and more challenging to apply correctly on continuous (rather than digital) signals.

And then, on top of the metrics mentioned above you have variations, like wavelet based cross-correlation instead of correlation, LZ complexity instead of Mutual Information and so on.

The field of Similarity Metrics is vast and there is a lot of work happening in the processing of Electroencephalography (EEG) signals in Functional Connectivity (which basically involves a lot of those evaluations of similarity under very challenging conditions of waveform noise and artifacts). For more information, please see this, this and this.

These publications are not exactly on your topic but they use a great deal of similar methods and tools for estimating the similarity between two time-series which is what you are trying to do as well. Therefore, when reading those, try to extrapolate the terminology and then try to search in remote sensing related journals for similar things. I am sure that such a search would return much more relevant results adapted to the specifics of estimating similarity of the laser pulse returns.

Hope this helps.