Integral over power spectral density The wikipedia entry on PSD has one confusing line:

  Summation or integration of the spectral components yields the total power (for a physical process) or variance (in a statistical process)

But would the total power correspond to the uncentralized second moment, so variance + mean2? In particular, there are stationary processes with a nonzero mean. Is the integral over the PSD still just the variance?

Summation or integration of the spectral components yields the total power (for a physical process)

yes, that's how "power density" is defined ...

  or variance (in a statistical process)

For any bound density, that's true.

  But would the total power correspond to the uncentralized second moment, so variance + mean²?

Yes, but if you have a process with a non-zero mean, that would imply that the spectral power density at $f=0$ diverges (a "Dirac delta $\delta(f)$" with "infinite height"), and that would imply the PSD isn't actually a function (it's still a distribution / test function, but seriously, this is just getting into math).

If you will, just imagine your "mean" is something at very low frequencies ($f\in[-\epsilon;\epsilon]$), but not concentrated to exactly $f=0$.

Still works out the same, for anything technical (because nothing will ever experience whatever process they're observing for eternity, so who's to say whether something happens at $0\,\text{Hz}$ or $10^{-24}\,\text{Hz}$?).

Honestly, the wording of the article isn't the greatest.


How to determine poles and zeros of the z-transform? This is a simple question, but I just don't understand how we determine the poles and zeros of a rational system function.

For example, for the LTI system described by this constant coefficient difference equation
$$
y[n]-\frac{5}{2}y[n-1]+y[n-2]=x[n]
$$
we can determine that

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

I understand why there are poles at $z=\frac{1}{2}$ and $z=2$, but I don't understand why there are two zeros at $z=0$. Even if I multiply $H(z)$ through by $z$, there would only be one zero, correct?

Just multiply numerator and denominator by $z^2$ to obtain

$$H(z)=\frac{z^2}{(z-\frac12)(z-2)}$$

from which you see that there's a double zero at $z=0$. Note that the number of zeros and poles is always equal if you include poles and zeros at infinity.


How to zoom an image using DCT transform I am familiar with zoom using fourrier transform on an Image.
I am also aware that DCT is used mainly to compress images.
Meanwhile, how can someone use this transform to zoom an image?

The following IEEE pdf discusses image interpolation (aka zooming) by DCT. Have a look at it
The main idea remains the same: the DCT spectra is zero padded and converted back into time to obtain the interpolated image which yields the zooming effect.


Find Amplitude of each period of time series I would like to determine the of each period amplitude of a time series (e.g. the amplitude of each day). The time series has a constant period of 1 day and varies only in amplitude. I have tried splitting the time series into individual periods then solving for the amplitude by re-arranging the equation:

$$y = A \cdot sin ( \frac{2 \pi}{\tau} \cdot t)$$

to:

$$A = \frac{y}{sin ( \frac{2 \pi}{\tau} \cdot t)}$$

where $A$ is amplitude, $\tau$. is the period and $t$ is time.

This works reasonably however there are issues with locating the beginning of each period.

Is there some more robust method I can use for this problem?

I assume from your solution that your time series closely fits a sine wave.  The method you are using is quite good if you restrict yourself to the set of data which is near the peaks.

The more standard way is to find a best fit sinusoid using linear algebra techniques.  You have stated the period, so you know the frequency.  Let's assume you don't know where the zero crossings are.  The problem then becomes to find the best values of $(a,b)$ so that

$$ Y = a C + b S $$

Where $Y$ is your signal values as a vector, $C$ is a cosine curve, and $S$ is a sine curve over one period.  Now dot this equation with $C$ and then $S$ to get:

$$ Y \cdot C = a C \cdot C + b S \cdot C $$
$$ Y \cdot S = a C \cdot S + b S \cdot S $$

The dot products are scalars.  Since $S$ and $C$ are orthogonal, their dot product is zero.  The dot product of $S$ or $C$ with itself is $N/2$ where $N$ is the sample count.  The equations then become

$$ Y \cdot C = a \cdot N/2 $$
$$ Y \cdot S = b \cdot N/2 $$

Solve for $a$ and $b$.

$$ a = \frac{ Y \cdot C }{ N/2 } $$
$$ b = \frac{ Y \cdot S }{ N/2 } $$

Your amplitude can now be found from $a$ and $b$.

$$ A = \sqrt{a^2+b^2} $$

This method is equivalent to a single bin of a Discrete Fourier Transform (DFT), so if you follow this then you understand how a DFT works.

To calculate the dot product, you take a summation:

$$ Y \cdot C = \sum_{n=0}^{N-1} { Y[n] C[n] } $$
$$ Y \cdot S = \sum_{n=0}^{N-1} { Y[n] S[n] } $$

That's all there is to it.


Stability condition for an iir filter I have been out of touch with DSP for a while. I was trying to refresh my basics and ended up with this doubt on IIR filter stability condition. I understood that for an iir filter to be stable, the poles should be inside the unit circle which was derived from its analog counterpart (i.e poles should be on Left hand of s plane). But then is there any extra condition that poles should be only real? Could it be complex also?

For a causal system to be stable, all the poles have to be inside the unit circle. In general, they may be complex, although they always come in complex conjugate pairs in order to keep the output real for real inputs.

So, the answer is no, there are no extra requirements for the poles to be on the real axis. However, they must be complex conjugates for the output to be real for a real input.

I would like to point out that the region of convergence of the Z Transform may include the unit circle, and yet the system may be non causal (but stable). For example, consider the impulse response: $h[n] = -a^{n}u[-n-1]$, for $|a| > 1$, which corresponds to a system having a pole outside the unit circle, and yet it is a stable (although non-causal) system.


detect and distinguish driving maneuvers via smartphone sensors I want to do some driving style analysis via smartphone.
Therefore, I programmed an app that records a 2D acceleration (when accelerating/braking and turning left/right). The app calibrates itself in a way that the initial misalignment between the reference systems of the car and the smartphone, that is fixed in the car while driving, is corrected.

For the start of my analysis, I want to focus on the longitudinal acceleration value. I want to distinguish between harsh acceleration and normal acceleration as well as harsh braking and normal braking (or even coasting). If I would be able to do that, I could determine over time whether a person is driving aggressive or not.

How can I distinguish harsh acceleration and normal acceleration as well as harsh braking and normal braking (or even coasting) in real-time?

I want to distinguish between harsh acceleration and normal acceleration as well as harsh braking and normal braking (or even coasting).

On the basis of accelerometer readings, you could place thresholds on acceleration regions or look for outliers based on a profile of "normal" driving.

A typical crash deceleration involves tens of gs (fig 7) and airbags deploy after -1.0g for more than 10ms (page 6).

Therefore, anything close to deviations of 1.0g is already looking suspicious in terms of "feelings" to be expected in normal driving conditions. (Of course, you have to make sure that the phone is fixed firmly on the car). As a side note, typical airliner flight does not seem to go too far above 1.0g (of course, it is more during take off and breaking at landing).

The other approach to take would be to carry your device along with you for a large number of "normal" trips and then process the collected waveforms to establish those "acceptable" limits yourself. This could be seen as Anomaly Detection. In other words, you would be deriving a set of features from your waveforms (modal mean of accelleration over a period of time, variance, median, typical waveforms involved with steering or breaking and accellerating, or other) and then evaluate how "close" a given driving session is to those metrics of normality.

  How can I distinguish harsh acceleration and normal acceleration as well as harsh braking and normal braking (or even coasting) in real-time?

Breaking and acceleration are covered above. By coasting, I am assuming that you mean allowing the car to roll in neutral. This is a rather difficult problem to be solved just by looking at a single accellerometer's output. It is very likely that it would have to be augmented by GPS or gyroscope sensors as well (most phones have all three types of sensors today).

To detect coasting, you would basically have to monitor very small changes in the attitude of the chassis. In constant velocity movement, the accelleration is zero and the engine is producing just enough power to counteract the friction of the tyres to the road and the chassis to the air. BUT! Part of the engine's torque is compressing (slightly) the suspension system. Therefore, a car under power tends to "lift its nose" (in side-view, the wheels want to turn in one direction, therefore the chassis tends to turn in the opposite direction. Try to notice the horizon in this video.). Of course, the road is not the track and for every torque setting, it is possible to find just the right downslope that would make the car appear to be coasting but it would still have a gear engaged. Therefore, inclination, either measured by gyro or GPS derived slope would be invaluable in determining coasting.

Hope this helps.


Remove artifact from space probe image This is a photo of our homeworld, taken by the Juno spacecraft recently, sling-shotting on its way to Jupiter. What it gained in speed, we lost in ours, however thankfully we will not be falling into the sun.

I think the South American continent is on the left.

However, we can notice that there is a sort of artifact across the image, a faint sort of blue bar that exists across the image. I am curious as to what may have caused that.

What I would really like to know though, is what image processing techniques us puny humans might muster in order to remove this artifact?

There is both a hue shift and a saturation shift from the full color image to the band artifact. Converting the color image to grayscale generates an image that still has an artifact, but that is less obtrusive.

There are "image inpainting" algorithms that you might try to correct the artifact, though it would help to first isolate the affected region. Here are a few thoughts on identifying the artifact:

Run Hough or RANSAC line fits for horizontal and vertical lines in the image. For a natural scene, and especially in a scene of the Earth, it would be odd to have horizontal or vertical lines corresponding to real features. Even the structures visible from space such as canals and large constructions tend to meander.
Try a form of image subtraction that subtracts a grayscale image from the original color image. Since the artifact is more grayscale-like and desaturated, you should be able to detect a region of hue and saturation differences. I would recommend working in HSV space for this operation. Run an algorithm in the vertical direction to look for a peak or valley in the resulting (color HSV - gray HSV) image.
Test row-by-row cumulative changes in Hue and Saturation. From row N to row N+1, the same X pixel generally won't change much in H, S, or V. However, you should notice a significant difference between the row just before the artifact and the row at the artifact. The difference might be as simple as the root mean square of H and S differences at each pixel, or even just the sum of differences. Or whatever. Pick a metric that is simple and that works.
Double check the original, raw image to determine whether the horizontal artifact extends into the black region of space on the left side of the image. In the copy of the image I downloaded, it appeared that the artifact did not extend into the black region.

Here are specs for JunoCam:
http://en.wikipedia.org/wiki/JunoCam

"The camera uses a Kodak image sensor, the KODAK KAI-2020, capable of color imaging at 1600 x 1200 pixels. It has a field of view of 18 x 3.4 degrees with three filters to provide color imaging."

That sensor should give a nice color image without weird artifacts when used on earth in gentler circumstances.

Although the artifact did not appear to extend into the black region of space at left, that could be a matter of dynamic range/windowing: there may be some difference in the charge on those pixels, but perhaps that difference is clipped to zero.

My guess is that this is a readout problem as the image rows are dumped. With older analog cameras we could see all sorts of crazy artifacts related to EM interference. With digital cameras that's usually much less of a concern, but one possibility is that some short, localized EM event occurred that somehow affected just certain rows.

If you zoom in to the top of the band artifact, you'll see desaturated rainbow artifacts. There are two non-adjacent rows that are brighter than their neighbors, too.

Since the artifact is lined up so nicely with the horizontal rows of the image, I would expect the problem is some internal electronics problem rather than a problem with lensing, color filters, flaring, or anything like that. Unless the camera became damaged by a particularly strong EM event, I would guess that the artifact would occur only rarely.


Detecting FSK near real time with c# I try to detect some FSK signals made by audio in near real time. I know that I will get some latency.
In the start I try to detect some tones of 10Khz in a duration of 75ms.
That is achieved by a Goertzel filter.
What is the prefered technique for that issue?
I thought of making a little buffer that will always be detected by the Goertzel, but this sounds bad because I don't know when the the signal starts.
Can you give me some help or idea ? Thanks

People have been decoding FSK signals for ages with very weak microprocessors.  The usual method is called zero crossing detection.

You basically just measure time (count the number of samples) between the times when the signal crosses the zero volt line.
The number of samples multiplied by the sampling rate tells you the rate as a time period.  One divided by that number is the frequency of your signal.

That works well enough with clean signals.  If your signal has noise mixed in with it, then you use a band pass filter to just pass the frequencies you are interested in.

I'd work up an FSK demodulator using zero crossing that works for clean signals.  Once that works, you can worry about filters and improving the performance.


Variance of function of random variable Is their an easier way to find variance of function of random variable?

Till now what I am doing is first find probability density function of (function of random variable) then integrate over range.

Assume $Y = g(X)$ be the function of RV $X$, then by using the following

$$E\{ g(X) \} = \int g(x) f_X(x) dx $$

variance of $Y$ can be computed without the computation of pdf $f_Y(y)$ as:

$$
\begin{align}
\text{Var(Y)} &= E\{ (Y-\mu_Y)^2 \} = E\{ Y^2 \} - (\mu_Y)^2 \\
& =  E\{ g^2(X) \} - E\{ g(X) \}^2 \\
& = \int g^2(x) f_x(x) dx  - \left(\int g(x) f_X(x) dx \right)^2 \\
\end{align}
$$


Can somebody explain SDF and MACE to me? I can't find any information on SDF and MACE save for a vague exert in Kinser's Image Operators books. Therefore I am dumbfounded. I can't understand what they do and what's their role in composite filters. Can somebody please explain them to me? Thanks.

Your question is allusive at that point: explanation of acronyms (SDF, MACE), need for composite filters... A starting point is  matched spatial filter (MSF). To overcome its limitations, several versions have been proposed. David Casasent is a common denominatoir of the ones you consider.

For the SDF filer, my reference is:

C. F. Hester and David Casasent, “Multivariant technique for multiclass pattern recognition,” Applied Optics, 1980

  A technique for multiclass optical pattern recognition of different
  perspective views of an object is described. Each multiclass
  representation of an object is described as an orthonormal basis
  function expansion, and a single averaged matched spatial filter is
  then produced from a weighted linear combination of these functions.
  The technique is demonstrated for a terminal missile guidance
  application using IR tank imagery.

with comment in:

J. B. Nelson, "Multivariant technique for multiclass pattern recognition: comment"

  In a recent article, Hester and Casasent discussed an interesting
  technique in multiclass pattern recognition. They mentioned that other
  approaches to the construction of an orthonormal basis function set
  were possible and specifically suggested the possibility of using a
  Karhunen-Loéve expansion as an alternative to their own Gram-Schmidt
  procedure. It is the purpose of this Letter to recommend yet another
  alternative, the singular value decomposition (SVD), which has very
  convenient conceptual and numerical advantages, and to show how it may
  serve as a foundation for future extensions of the theory. The scope
  of the present discussion is restricted solely to the digital
  processing aspects of Ref. 1. It will begin with a brief review of the
  problem in the slightly modified matrix notation that is more
  compatible with the SVD technique that follows.

Abhijit Mahalanobis, B. V. K. Vijaya Kumar, and David Casasent, Minimum average correlation energy filters, Applied Optics, 1987

  The synthesis of a new category of spatial filters that produces sharp
  output correlation peaks with controlled peak values is considered.
  The sharp nature of the correlation peak is the major feature
  emphasized, since it facilitates target detection. Since these filters
  minimize the average correlation plane energy as the first step in
  filter synthesis, we refer to them as minimum average correlation
  energy filters. Experimental laboratory results from optical
  implementation of the filters are also presented and discussed


Does a matched filter really improve SNR? Does the matched filter actually improve SNR or does it maximize it (making sure that it is not less than what it should be) ? improving SNR would not make sense to me as suppose a receiver has a certain NF that defines an SNR for a signal at a certain power level and a certain BW, improving SNR would mean reducing NF, and that wouldn't make sense, would it?

Assuming a pulse shaped  BPSK signal with a symbol or bit rate =1/BW (a roll off almost equal to zero), the SNR at the matched filter output =Eb/No, multiplying the denominator and the numerator by BW we get: signal power/ noise power and noise power would of course be a function of NF, so in that case where is the improvement? or do you want to say that already assuming that signal BW=1/rate because of pulse shaping  is the condition that did the job?

There is no contradiction. Matched filters are often used because they are an implementation of a sliding correlator, which is the optimum detector (in the maximum-likelihood sense) for linear modulation in the AWGN channel. That is, the SNR at the output of the matched filter will be maximized when the filter lines up with the point of reception of the signal that the filter matches.

This doesn't contradict the known properties of your receiver (i.e. its noise figure). The point to notice is that the SNR can be different at the input of the matched filter than it is at the filter output. The SNR at the matched filter input (i.e. your front end's output) is determined by the signal power at the receiver input and the noise figure of the front end. But by applying signal processing techniques, like matched filtering, you can improve that SNR to achieve more accurate detection.


what happens If I change the phase and magnitude of a signal I have a speech signal. I want to know what happens when I make some changes on the phase of its frequency reaponse. When I set its phase to zero I observed that the new signal is even and the amplitude of the first sample became larger. Does it mean all the frequency components get the zero delay? But there is again amplitude in middle part of the signal. Is there any better interpretation for makeing the phase of signal zero?

When I set the frequency resposnse to the constant value of its average of magnitude what will extactly happen?

Thanks in advance for any help

Regarding the phase: if you set all frequency components to have zero phase, then the signal becomes essentially a sum of sinusoids that tends to an impulse as the frequency increases. In other words,

$$\delta(t) = \int_{-\infty}^\infty \cos(2\pi ft) df.$$

In your case, the maximum frequency is limited, which means that the signal will actually become a sinc. Try this in Matlab (or similar):

fs = 1000;
t = -0.1:1/fs:0.1;
f = 0:200;

s = zeros(1,length(t));
for ff = f
    s = s + cos(2*pi*ff*t);
end
plot(s);

and you'll get a plot like this:

When each cosine has a different amplitude, the signal still will tend to acquire a sinc-like shape.


Aliasing and DTFT of a real signal We are analyzing a real signal with the DTFT. Since we are using a limited number of samples it's like we are transforming a finite signal.

As I remember, the FT of a finite signal has an infinite band spectrum. So I suspect that the DTFT of a sampled signal (with a finite number of samples) always produce aliasing (since in the DTFT spectrum each replica is an infinite band). Is this correct?

We are analyzing a real signal with the DTFT. Since we are using a limited number of samples it's like we are transforming a finite signal.

No. The DTFT takes an infinite discrete time signal as an input but the spectrum is continuous. I think you are confusing the DTFT and the DFT. They are different things.

  As I remember, the FT of a finite signal has an infinite band spectrum.

Yes

  So I suspect that the DTFT of a sampled signal (with a finite number of samples) always produce aliasing (since in the DTFT spectrum each replica is an infinite band). Is this correct?

Yes and no. Sampling a finite signal does with a finite sampling frequency will always induce some amount of aliasing, although in practice you can make it almost always "small enough" so it's acceptable for your application. That has nothing to do with DTFT (or DFT). Aliasing is caused by sampling, not during the analysis of the sampled data. IN other words, the aliasing happens regardless of whether you apply a DTFT or not.


Optimized resource allocation problem I am from ECE background and trying to solve channel allocation problem. Let's assume I have three users and three available channels. I would like to allocate channel among them in such a way that average BER would be minimum.
For example, consider BER matrix for each user for each frequency allocation
\begin{bmatrix} User/frequency & User1 & User2 &  User3 \\ f1 &0.01 & 0.05 & 0.3\\f2 & 0.12 & 0.045 & 0.02 \\ f3 &0.02 & 0.01 & 0.01 \end{bmatrix}

Now, I would like to allocate channels f1--> User1, f2--> User3 and f3-->User2, which gives least average BER (Sum of each users BER/no of users). (Note: Same frequency can't allocated to different users.)

Brute force is obliviously one solution, but I am trying to find optimized solution.

Your problem is called an assignment problem. Check out the Hungarian algorithm. Instead of average BER cost you can minimize total cost as it gives the same optimal solution.


extend parts in an existing audio file to make it sound longer without stretching the file I'm looking for a way to extend certain parts of an exiting audio file to make these parts sound like the where played longer. And I dont want to change play speed or stretch the file.
Lets say u have a digital sound file of someone playing a trumpet. And he's playing diffrent tones. Now I want to jump to a certain point in this file and hold this very tone as long as I want. To make it sound like it was played very very long. Without anyone noticing that it was artificialy extended.
And I want to do this with not just recordings of instruments.

What technics do you suggest?

Maybe something like granular synthesis or any other synthesis technic?

It would be so awesome if someone could help me with this.
Thanks in advance

Sounds for me, like non-linear time stretching can accomplish what you want. In your example, you would apply a stretching function with a slope of one most of the time and for the single note (you want to be played very long) the slope then needs to be strongly different from one.

Overview paper: Jonathan Driedger and Meinard Müller. A Review on Time-Scale Modification of Music Signals, Applied Sciences, 6(2): 57–82, 2016. [DOI]

Code: https://www.audiolabs-erlangen.de/resources/MIR/TSMtoolbox


How to derive 2nd order Butterworth condition for the damping coefficient mathematically? The magnitude of 2nd order low pass filter is given as
$$|H(\omega)|^2= \frac{1}{(1-(\frac{\omega}{\omega_o})^2)^2+(\frac{2\zeta\omega}{\omega_o})^2}$$
Now in order to achieve maximally flat within pass band, we take the derivative of this equation and set it to zero to find the extremums.
\begin{equation*}
    \frac{\partial(\frac{1}{|H(\omega)|^2})}{\partial\omega} = \frac{\partial[(1-\omega^2)^2+(2\zeta\omega)^2]}{\partial\omega}
\end{equation*}
\begin{equation*}
\begin{split}
    2(1-\omega^2)(-2\omega)+8 \zeta^2\omega &= 0 \\
    \omega(\omega^2-1)+2\zeta^2\omega &= 0 \\
    \omega (\omega^2-1+2\zeta^2) &=0 \\
    \end{split}
\end{equation*}
Then we get $\omega =0$ or $\omega^2 = 1-2\zeta^2$
Now my question is, why do we set the second root to be zero to obtain no ripples in the passband so that we can derive:
$$\zeta=1/\sqrt{2}$$
I understand that the slope of the magnitude is zero at these roots, but I couldn't really interpret the idea here. Any help is appreciated.

You have to check the next derivative, because the first one is zero anyway at $\omega=0$. The second derivative is
$$3\omega^2+2\zeta^2-1\Big|_{\omega=0}=2\zeta^2-1\tag{1}$$
which can be made zero by choosing $\zeta=1/\sqrt{2}$.


Fourier Transform of exponential While solving Example 4.1 of Signals and Systems by Alan Oppenheim. Example 4.1 is:
$$ x(t)=e^{-at}u(t), a>0$$
and the transform I get is:
$$ X(j\omega)\frac{1}{a+j\omega}, a>0$$

The problem is understanding the sketch:
The sketch is drawn for both the values $+a$ and $-a$> i dont understand why it $-a$ when $a>0$ is defined.

The plot is of $$\mid X\left(i\omega\right) \mid = \sqrt{\left(\frac{1}{a+j\omega}\right)\left(\frac{1}{a-j\omega}\right)} = \frac{1}{\sqrt{a^2 + \omega^2}}$$

against $\omega$

In particular $\omega$ can be equal to $-a$.

This checks out with Wolfram alpha


Understanding legacy code: Algorithm to remove radial lens distortion The method below to remove lens distortion from a camera was written more than ten years ago and I am trying to understand how the approximation works.

void Distortion::removeLensDistortion(const V2 &dist, V2 &ideal) const
{
 const int ITERATIONS=10;
 V2 xn;
 xn.x=(dist.x - lensParam.centerX)/lensParam.fX;
 xn.y=(dist.y - lensParam.centerY)/lensParam.fY;

 V2 x=xn;
 for (int i=0;i<ITERATIONS;i++) {
   double r2 = Utilities::square(x.x)+Utilities::square(x.y);
   double r4 = Utilities::square(r2);
   double rad=1+lensParam.kc1 * r2 + lensParam.kc2 * r4;
   x.x/=rad;
   x.y/=rad;
 }
 ideal.x=x.x*lensParam.fX+lensParam.centerX;
 ideal.y=x.y*lensParam.fY+lensParam.centerY;
}

As a reminder:

lensParam.centerX and lensParam.centerY is the principal point
lensParam.fX and lensParam.fY is the focal length in pixel
lensParam.kc1 and lensParam.kc2 are the first two radial distortion coefficients. This is $k_1$ and $k_2$ in the formula below.

The formula to add lens distortion given the first two radial distortion parameters is as follows:

$$
x_{distorted} = x_{undistorted} \cdot (1+k_1\cdot r^2 + k_2 \cdot r^4) \\
y_{distorted} = y_{undistorted} \cdot (1+k_1\cdot r^2 + k_2 \cdot r^4)
$$
where $r^2=x_{undistorted}^2+y_{undistorted}^2$ and $r^4=(r^2)^2$

In the code above, the term $(1+k_1\cdot r^2 + k_2 \cdot r^4)$ is calculated and saved in the variable rad and the distorted x is divided by rad in each of the ten iterations.

All of the cameras we use have a pincushion distortion (so $k_1<0$)

The question is how is this algorithm approximating the undistorted image points?
Do you know if there is any paper in which this algorithm is proposed?

The opencv undistortion may be a bit similar, so that link may be useful but it is not quite the same though.

Welcome to SP. I happen to have already answered your question in SO: see point 2 of my answer here.


How can i measure shape elongation feature in matlab or python opencv? How can i measure shape elongation faeture in matlab or python opencv?
There is a function library or i have to implement to my own?
In case there is no function, i found this definition:

"The less-common elongation shape factor is defined as the square root of the ratio of the two second moments in of the particle around its principal axes"

f_elong = sqrt (i2/i1)

Where i found these first and second moments?

In the Image Processing Toolbox for MATLAB there is a function called regionprops, which may be what you need. It takes a binary image, finds connected components in it, and returns various measurements of each connected component.  Here's a list of things it can compute:

Shape Measurements

  'Area'              'EulerNumber'       'Orientation'
  'BoundingBox'       'Extent'            'Perimeter'
  'Centroid'          'Extrema'           'PixelIdxList'
  'ConvexArea'        'FilledArea'        'PixelList'
  'ConvexHull'        'FilledImage'       'Solidity'
  'ConvexImage'       'Image'             'SubarrayIdx'
  'Eccentricity'      'MajorAxisLength'
  'EquivDiameter'     'MinorAxisLength'

Pixel Value Measurements (requires grayscale image as an input)

  'MaxIntensity'
  'MeanIntensity'
  'MinIntensity'
  'PixelValues'
  'WeightedCentroid'


Size (bytes) of a sampled signal I am a beginner and maybe I have the false idea but I will thank your understanding, my background is not so on the DSP side.

I have a KX134-1211 accelerometer that delivers a digital output.
It is triaxial 16 bit sensor and delivers at 2 axis 8500 Hz as signal bandwidth.

This means the sensor fmax is 8500 Hz right? Since it is already digital how much sampling rate do I need to get the maximum signal banwidth, 8500 or 17000 Hz?

If I sameple the sensor for at 8500 during 1 second I will get 8500 values, how much bytes is this, I want to calculate how much size each recording will take.

Please have a look at Page 5 of the DataSheet here.

It says the signal Bandwidth is 8500Hz only when output data rate is 25600Hz and Low Pass Filter setting = 1. And that the Bandwidth changes based on Output Data Rate.

So, I think you do not have to worry about sampling rate, since ODR(Output Data Rate) is 25600Hz and that means 16-bits(2Bytes) readings coming at this rate.

So, 1 second of recording will require : 1*25600*2Bytes = 51200/1024 KB = 50KB, of memory space. And, I am not sure but I think this is just for X-Axis measurements by the accelerometer.


Estimate smartphone accelerometer bias I was planning to develop Android / iOS applications that enable users to measure 3D length using their smartphones.

According to this question, you need to know at least the time-varying bias that affects the accelerometer's signal. Hence, I need to estimate the bias every time a user is going to measure. How can I do that? I've seen articles about this estimation but to be honest I'm not expert in signal processing so I didn't understand them. Can you explain a method to estimate?

According to this IEEE article: https://github.com/liulinbo/slam/blob/master/A%20robust%20and%20easy%20to%20implement%20method%20for%20IMU%20calibration%20without%20external%20equipments.pdf

You can model errors this way:
$$ a = f*a' + g + b + \eta$$
where $a$ is the actual acceleration, $f$ is a 3x3 matrix to model scaling, misalignments, cross-axis and ... errors. $a'$ is sensor's data, $g$ is gravity, $b$ is 3x1 matrix to model bias, and $\eta$ is 3x1 matrix to model noise.

You can calibrate accelerometer by reading $a'$ from sensor data when the accelerometer is in static positions. In static positions, the only force effecting accelerometer is gravity. So by minimizing this summation, $f$ and $b$ can be computed(calibration):
$$ (|a_0|^2 - |g|^2)^2 + (|a_1|^2 - |g|^2)^2 + ... + (|a_N|^2 - |g|^2)^2 $$
where $N$ is number of static positions.


How to Apply Window Functions in Frequency Domain? I am trying to apply a window function to smooth my signal and finally obtain the peaks of the signal. Yet, for spectral leakage I apply window function in time domain. After taking the Fourier amplitude spectrum ı want to further smooth the signal. Yet I do not know how to apply smoothing windows for frequency domain. Should I use convolution for the Fourier amplitude spectrum of the window and Fourier amplitude spectrum of the data?

What is the meaning of smoothing window length with 0.1 Hz or 0.2 Hz?

You need to use the Convolution Theorem.

Namely, when you want to apply an LPF Kernel in one domain (Using convolution) you can use multiplication in the other domain.

The tricky part is the dimensions of the signals (Since we're dealing with Discrete signals).
Since the LPF Kernel in the convolution domain might be different then the size of the signal it is applied to, so how do we multiply?

Well, the solution is given by interpolation to the correct dimension.
This is done by adding zeros at the end of the signals.
Pay attention that the end of the signal in the frequency domain is tricky.
The end are the samples before the sampling frequency (On both sides as it is symmetric), add zeros and keep the signal symmetric (If you signals are real in the other domain).


Difference of doing a PSD estimate of data and logarithmic transformed data? What is the difference between doing a PSD estimate with data and the same data but which is logarithmically transformed before the estimate? Does it make the data more sinusoidal in nature?

An exercise in a book asks this question:

"For the lynx data, compare your spectral analysis results from the original data, and the data transformed first by taking the logarithm of each sample and then by subtracting the sample mean of this logarithmic data. Does the logarithmic transformation make the data more sinusoidal in nature?"

The exercise is instructive. Since this is an exercise, I will not do the full problem.  I leave that to you.

grabbing the lynx data from

https://www.encyclopediaofmath.org/index.php/Canadian_lynx_data.

The data is counts and there are no years with a zero count.

and putting it into matlab and taking the log and subtracting the mean of the log,

Would you say that the data looks sine like ?

This technique works on certain kinds of data.  Think of it as a tool and is covered in a number of basic stats book.


What is the effect of noise on Shannon entropy? Since entropy is a measure of uncertainty or randomness, intuitively we would suppose that adding noise to an image would increase its entropy since we are now more uncertain about the information of the image. I actually confirmed (partially) this assumption with the following MATLAB code where i noticed that when adding gaussian noise to a grayscale image, the entropy increases, whereas adding salt & pepper noise slightly decreases the entropy. Why is this happening and how does noise affect overall the Shannon entropy of an image?

This is my MATLAB code:

load('trees.mat');
RGB=ind2rgb(X,map);
GRAY=rgb2gray(RGB);
vecRGB=RGB(:);                                 %--vector of RGB image
vecGRAY=GRAY(:);                               %--vector of GRAYSCALE image
HVECR=entropy(vecRGB);                         %--entropy of RGB image
HVECG=entropy(vecGRAY);                        %--entropy of GRAYSCALE image
GRAY_GAUSSIAN=imnoise(GRAY,'gaussian');        %--add gaussian noise to GRAYSCALE image
GRAY_SALT=imnoise(GRAY,'salt & pepper');       %--add salt & pepper noise
vec_GRAYG=GRAY_GAUSSIAN(:);                    %--vectorize image with gaussian noise
vec_GRAYSP=GRAY_SALT(:);                       %--vectorize image with salt & pepper noise
HVECGG=entropy(vec_GRAYG);                     %--entropy with gaussian noise
HVECGSP=entropy(vec_GRAYSP);                   %--entropy with salt & pepper noise

And these are the entropies i got:

Entropy of Grayscale Image (no noise)=5.4723
Entropy of Grayscale    Image (Gaussian Noise)=7.6948
Entropy of Grayscale Image (Salt &    Pepper Noise)=5.4380

The grayscale image is actually a discrete image. entropy calculates a histogram and from that extracts "empiric" probabilities, to be used in the common information formula

$$I= -\log_2 (P(g))$$

with $P(g)$ being the (empiric) probability of an intensity value $g$.

Now, entropy $H$ is the expectation of Information. For your whole image $V$ that means

$$\begin{align}
H(V) &= E\left\{ I(g_{x,y})\right\}\\
&= \sum\limits_{x,y\in \text{idx(V)}} P(g_{x,y})\cdot\left(-\log_2 (P(g_{x,y}))\right)\\
&= -\sum\limits_{x,y\in \text{idx(V)}} P(g_{x,y})\log_2 (P(g_{x,y}))
\end{align}$$

If you look at the product in the sum intensely enough, you'll notice that the sum will take a maximum value if all properties are equal – that is, if every intensity value is just as likely as every other. The histogram of that image would be a flat line.

Adding "well-behaved" noise will make almost any image look more like that – in the end, pixels with identical intensity will be relatively common in most natural images, and adding noise will "break" that.

Now, salt'n'pepper is different: you add extreme values to your image. That will make the intensity at many pixels clip at the maximum intensity (or the minimum one), skewing the histogram away from flat to a very "horned" thing, which has less entropy.


Design of a low pass filter with Red Pitaya's FPGA I'm trying to design a low pass filter in Verilog for use on a Red Pitaya's FPGA (Xilinx® Zynq®-7010).
All I'm trying to get the low pass filter to do is to obtain the DC component in the signal.
The only part of the signal I require is the DC component.
Can someone put me through on how to go about the design as I am new with Red Pitaya?

There are 2 steps for the solution
Step 1 - Design some kind of low-pass filter to measure the DC component.
In order to perform step 1 you need to answer a few questions :
What do you mean by DC? For example, measure the DC for the last 10 seconds?
What is the input sampling rate? What is the desired output sampling rate?
You need to provide this information, we cannot guess it for you.
Step 2 - Implement the filter in Verilog.
There are examples on opencores. I can also give you a hand if you need help with implementing a FIR or IIR filter in Verilog.
You could also look at the Xilinx FIR Cores in the Core generator.


Frequency translation of an ideal low pass filter I am trying to create a discrete-time filter with a Fourier transform as follows.

$X(\omega) =
\begin{cases}
1, & T - W \leq \omega \leq T + W\\
0, & \text{all other values of } \omega
\end{cases}, $

Looking at the provided solution, I have understood that the idea is to create a low-pass filter and then translate it in the frequency domain by $T$.

For the low-pass filter I've got the following time-domain function.

$x(n) = \frac{1}{\pi}\frac{\sin(Wn)}{n}$

To translate a function in frequency, it can be multiplied in the time domain by $e^{jTn}$. In the solution, they have first stated that this is correct, but in the final answer they have instead multiplied it by $2\cos(Tn)$. They also wrote that $2\cos(Tn) = e^{jTn}+e^{-jTn}$, which I have no problem with, but I don't understand why they went for this instead of $e^{jTn}$.

So to be clear, my question is:  why did they multiply the low-pass filter by $2\cos(Tn)$ instead of $e^{jTn}$ to translate it in the frequency domain?

Multiplying the impulse response by $e^{jTn}$ would give a complex impulse response and a complex-valued filter output. Seems that they were aiming for a real-valued output. They did this by summing two filters, one with the pass band at $T - W \le \omega \le T + W$ and the other with the pass band at $-T - W \le \omega \le -T + W$. If there is overlap between the two pass bands then the frequency response of the composite filter will have a weird shape like:

Seems that they did not intend there to be overlap, but that $T \gt W$.


Relation between sawtooth Fourier coefficients and its DFT I'm having some trouble with understanding the DFT of a sawtooth single period signal and its relation with sawtooth Fourier coefficients.

Let's say I have a signal
$$
s(t) = \frac{At}{T} - \frac{A}{2} \qquad t\in[0,T)
$$

If I plot its DFT real and imaginary components I have something like this

Where $A=10$, $T=1$, $N=50$ (number of samples) and the DFT is normalized by $\tau = T/N$ and $1/N/\tau$ for the inverse.

Now, the imaginary part as far as I can understand comes from the Fourier coefficients
$$
c_n = \frac{iA}{2\pi N}
$$
These coefficients (see e.g. Brigham, section 5.2) equal the fourier transform $S(t)$ scaled by $\frac{1}{T}$ and sampled at $n/T$. Given we have a discrete signal, we need to account for aliasing in the frequency space, which can be represented by repeating $c_n$ at $N$ intervals and summing with themselves, something like the following? probably missing a constant
$$
Im(S_n) = \sum_{m=-\infty}^{\infty} c_n\delta(n-mN) \qquad n \in [0,N-1]
$$

What about the non zero real part? As an afterthought I can rationalize it because the sawtooth is non zero at $t=0$, so it must have a non zero mean DFT. But if I take a sawtooth centered in $-T/2,T/2$ it's still there even if the signal crosses the origin.

I'd like to find a better relation that gets me from the Fourier series to the DFT, accounting both for the imaginary and the real part. Any pointer to good resources to help understand the whole thing would be more than welcome.

Just looking at the DFT, the real part, with my normalization should be:
$$
-\frac{A}{2}\tau\sum_{n=-\infty}^{\infty}\delta\left(f-\frac{n}{T}\right)
$$
How do I justify this result? especially the multiplication factor?

(sorry if the math is a bit flaky, I'm still trying to understand everything about DFT, Dirac combs and most of all normalization constants, every single source has them different...)

Let me use the definitions of the Fourier series and the DFT as they are commonly used in signal processing. The Fourier series of a $T$-periodic function $f(t)$ is given by

$$f(t)=\sum_{k=-\infty}^{\infty}c_ke^{j2\pi kt/T}\tag{1}$$

where $c_k$ are the complex Fourier coefficients:

$$c_k=\frac{1}{T}\int_0^Tf(t)e^{-j2\pi kt/T}dt\tag{2}$$

The length-$N$ DFT of a sequence $x[n]$ is defined by

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}\tag{3}$$

and the inverse DFT is

$$x[n]=\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j2\pi nk/N}\tag{4}$$

If you sample $f(t)$ at $t=n\tau$ with $\tau=T/N$ you can approximate the integral $(2)$ by the following sum:

$$c_k\approx\frac{1}{T}\sum_{n=0}^{N-1}f(n\tau)e^{-j2\pi kn\tau /T}\tau\tag{5}$$

With $x[n]=f(n\tau)$ and $\tau=T/N$, Eq. $(5)$ can be rewritten as

$$c_k\approx\frac{1}{N}\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}=\frac{1}{N}X[k],\qquad |k|\le\left\lfloor \frac{N}{2}\right\rfloor\tag{6}$$

Equation $(6)$ says that the Fourier series coefficients $c_k$ of the $T$-periodic function $f(t)$ are approximated by the (scaled) DFT coefficients of the sequence $x[n]=f(n\tau)$, with $\tau=T/N$.

Of course this is also the case for your function. With

$$f(t)=A\left(\frac{t}{T}-\frac12\right),\qquad 0\le t< T\tag{7} $$

and $f(t)=f(t+T)$, we get

$$x[n]=f(nT/N)=A\left(\frac{n}{N}-\frac12\right),\qquad 0\le n< N \tag{8}$$

The Fourier coefficients of $f(t)$ are

$$c_k=\frac{jA}{2\pi k}\tag{9}$$

The DFT of $x[n]$ is given by

$$X[k]=\begin{cases}-\frac{A}{2},&\quad k=0\\
\frac{A}{e^{-j2\pi k/N}-1} = e^{j\pi k/N} \frac{jA/2}{\sin(\pi k/N)}=\frac{A}{2}(-1+j\cot(\pi k/N)),&\quad 0<k<N
\end{cases}\tag{10}$$

The first thing to notice is that $f(t)$ is an odd function (with a discontinuity at $t=0$), and, consequently, its Fourier coefficients are purely imaginary. On the other other hand, the sampled version of $f(t)$ is not odd because $x[0]\neq 0$; that's why we get a (relatively small) non-zero real part in the DFT coefficients. However, as predicted by Eq. $(6)$, the scaled DFT coefficients $X[k]/N$ still approximate the Fourier coefficients quite well. The zero real part of $c_k$ is "approximated" by the value $-A/2N$, which is small for large $N$. The figure below shows the imaginary part of $X[k]/N$ (red) and the (imaginary part of the) Fourier coefficients $c_k$ (blue):

The increasing difference between $c_k$ and $X[k]/N$ with increasing $k$ is due to aliasing caused by sampling $f(t)$. Note that $f(t)$ is not band-limited, so you will always have some amount of aliasing, no matter how fast you sample. The difference in the real part ($-A/2N$ vs $0$) is caused by the fact that $f(t)$ is odd, but $x[n]$ isn't because $x[0]=f(0^+)=A/2\neq 0$.

EDIT: Here comes the derivation of the DFT of $x[n]$ defined by $(8)$. From $(3)$ we have

$$X[k]=\frac{A}{N}\sum_{n=0}^{N-1}ne^{-j2\pi nk/N}-\frac{A}{2}\sum_{n=0}^{N-1}e^{-j2\pi nk/N}\tag{a}$$

Let's first consider the case $k=0$:

$$\begin{align}X[0]&=\frac{A}{N}\sum_{n=0}^{N-1}n-\frac{A}{2}\sum_{n=0}^{N-1}1\\&=
\frac{A}{N}\frac{(N-1)N}{2}-\frac{AN}{2}\\&=-\frac{A}{2}\end{align}\tag{b}$$

For $k\neq 0$, the second term in Eq. (a) equals zero. From the formula

$$\sum_{n=0}^{N-1}nq^n=\frac{(N-1)q^{N+1}-Nq^N+q}{(1-q)^2},\qquad q\neq 1\tag{c}$$

(which can be obtained by taking the derivative of the standard formula for a geometric series), and by noting that $e^{-j2\pi k}=1$, we get after some manipulations

$$\begin{align}X[k]&=\frac{A}{N}\sum_{n=0}^{N-1}ne^{-j2\pi nk/N}\\&=\ldots\\&=\frac{A}{e^{-j2\pi k/N}-1},\qquad k\neq 0\end{align}\tag{d}$$

from which all other representations shown in Eq. $(10)$ can be derived quite easily.


Fully random PN sequence Can a CDMA PN sequence be fully random, as in encryption-grade randomness, or must it be deterministic? I have tried the first bytes of SHA-256 hashes as PN sequences, but somehow they don't coexist too well.

I know a deterministic code can be generated on-the-fly so no storage is needed, but let's say I have 500 fully random PN sequences stored on a device. Is that good or do you really need deterministic codes like the Gold codes?

What do you mean under encryption-grade randomness?

It is not recommended to use linear random generator in most of the encryption cases.

PN-Sequence is a mutation of the LFSR. Because of the linearity, all future values can be calculated easily. So be careful using linear pseudo random generators.


Zero-phase vs Linear-phase vs Nonlinear-phase System I don't understand what is the difference between these. I'm trying to figure out how to classify this system's phase response:

Would someone please tell me what it is and how you can tell?

The phase seems to be either 0 or $\pi$, which correspond to real positive and real negative frequency responses respectively.

So this can be considered a real (zero phase) system.


Compressed Sensing and Sparsity Ok for first information I would like to let you know that I originally do not have any formal exposure in signal processing so I would appreciate if you consider me as a beginner. In particular in this project I have to deal with the method of compressed sensing (CS). As some of you might know that CS relies on the sparsity of the signal in some basis.

My first question is, when a reference mentions "wavelet" with no specification of which type of wavelet being referred to, what would be the type of wavelet which is common enough so that this reference doesn't mention it, if there is any?

Second, the main matrix equation in CS is
$$
\begin{align}
y &= \Phi \Psi s\\
\Updownarrow\\
y &= \Phi x,\quad x = \Psi s
\end{align}
$$
where $y$ the measured samples, $\Phi$ measurement matrix, $\Psi$ the representing basis matrix, and $s$ the expansion coefficients in $\{\Psi_i\}$ basis which we want to solve for. Ok up to this point I hope I haven't made any mistake in interpreting that equation.

My question is how do people normally chose the basis? I think this choice of basis is critical because the expansion coefficients $s$ depends on the chosen basis.
As we would want $s$ to be as sparse as possible, doesn't it also mean that we need to carefully choose the most favorable basis?
If so, how would I decide the right choice? I read in some references they mostly use wavelets representation as the basis.

and by the way, the mathematical expression of the expansion in wavelets domain contains all possible indices of the bases (wavelets), that is the signal is the sum over both indices (dilation and translation indices) of the bases.

How should I choose which pairs of indices I better put in the matrix $\Psi$, if for example the signal is $n$ elements long? In that case I obviously only need $n$ bases and the corresponding coefficients out of the infinite possibility.
Can I randomly take any $n$ bases provided any two different bases are always orthogonal?

As far as choosing the right dictionary (basis) goes, there are several things to take into account. You do want your representation $s$ to be as sparse as possible and that depends a lot on the chosen dictionary. In practice, people are often looking at a particular type of signals $x$ which they know from domain knowledge to be sparse in some dictionary (discrete wavelet transform, discrete Fourier transform, etc.). However, you also need the dictionary $\Psi$ to "work well" together with the measurement matrix $\Phi$. This can be measured by their coherence (see e.g. Candès & Wakin, 2008). If you are free to choose $\Phi$ as you like, you can select the dictionary $\Psi$ that gives you the sparsest representation and then select $\Phi$ to be incoherent with $\Psi$. In practice $\Phi$ is often limited by hardware constraints and then you may have to select a suboptimal $\Psi$ with respect to sparsity in order to achieve incoherence with $\Phi$.

Another approach is dictionary learning where you can use (many) examples of the type of signal $x$ you are looking for to learn a dictionary $\Psi$, see for example Aharon, Elad, Bruckstein, 2006.


Transfer function between voltage and decibel output for in-ear headphones What would be a decent formula for approximating the connection between voltage and decibel output for a typical set of unpowered in-ear headphones?

I'm generating sine waves from the computer: an array of floats between -1.0 and +1.0 gets sent out through the speaker:

X_out = h * sin( f * 2π * t ); // h is the amplitude

Let's say we assume the user to ramp up their volume so that h=1 corresponds to a loud  but comfortable listening level.  So:

h=1 ~ dB_out=80dB

Now, how can I mathematically connect h and dB_out?

I've read that:

v_out ~ h

Somewhere I have seen the formula:

SPL = 20 * Log10( v_out / v_ref )

So, rearranging:

v_out = v_ref * 10 ^ ( SPL / 20 )
k * h = v_ref * 10 ^ ( SPL / 20 )
h = k' * 10 ^ (SPL / 20 )
  = A * exp( SPL ) + B,  for suitable A, B

Is this working valid?

Here is the alternative formulation of the question:

Given a frequency f, and a desired decibel level,  what h do I need to generate that decibel level?

That depends to a very large degree on the type & sensitivity on the headset, the gain of the amplifier and the frequency. I don't think there is a reasonable rule of thumb, this will be all over the place.


Constraint on a filter First of all thanks for your patience: it is the first time for me posting a question in this forum. I am not a DSP expert, but I should get by if you give me in depth explanation.

This is an example of my raw data

and a link to the actual dataset (time=t0, acceleration=acc).

My problem is the following: I am measuring accelerometer data for an experiment, but I am having a hard time filtering out the raw acceleration data for double integration of velocity and position. The first thing that I do is to remove the dc component of the data via

#python code
func = interp1d(t0,acc)
dt = 10./1000. # [secs]
ta = np.arange(t0[0],t0[-1],dt)
a = func(ta)
n = len(a)
freq = np.fft.fftfreq(n, d=dt)
freq[0]=1e-16
F = np.fft.fft(a)
F[0]=0
new_acc = np.real(np.fft.ifft(F))

In particular, I have the physical constraint the integrated velocity should be always positive, or equal to zero where the acceleration is equal to zero.
However, when I integrate the new_acc,

v = cumtrapz(new_acc)*dt

this is what I get:

Leaving alone for a minute the position (which is ultimately really what I am interested in)
my question is:

is there any way to design a filter such that the integral of the filtered signal is always >=0 ?

Alternatively, if you believe that the one above is the wrong question to ask:

what is the best way to integrate (twice) this dataset?

Beware.  A Gaussian random walk has a non-zero distribution width that increases with time.  Thus the double integral of that random walk will likely very quickly zoom off the edge of your solution space.  A linear filtering of that random walk won't behave any better.  Whether or not mixed with the "real" acceleration signal.


Why can convolution only be applied to compute the output of a linear filter? We apparently cannot compute the output of a bilateral filter (BF) using convolution (with the image) because the BF is a non-linear filter. In general, why can convolution only be applied to compute the output of a linear filter?

This is a good question and something that I remember asking myself when I first learned about impulse responses and convolution.

To understand this, it is first necessary to understand the significance of impulses and impulse responses. Referring to the image below, you can see that an impulse is an instantaneous like input and the impulse response is the decaying output.

So why is the impulse input significant? It's significant because you represent any arbitrary signal input as an array of impulses! Literally any input you could ever want to feed into a filter can be thought of as an array of impulse inputs.

Now that we understand the significance of impulse input, what is the significance of the impulse response? Well for linear systems we know that the output must follow certain properties for given inputs:

Additivity: y(a + b) = y(a) + y(b)
Homogeneity: y(cx) = c y(x)

What this tells us is that if we feed in multiple inputs (i.e. feed in an array of impulse inputs), we should be able to easily compute the output as a summation time shifted impulses responses. Convolution is simply the mathematical operation that performs this process. If you look at the diagram below, you can see that the convolution operation is simply summing together time shifted impulse responses.

So to answer your question as to why convolution only works on linear filters. It comes down to the fact that convolution relies on the two linear properties (additivity and homogeneity) in order to be accurate in predicting an output. If these two properties aren't true for a system, then the impulse responses can not be summed to calculate the output thus "breaking" the usefulness of convolution.

Although my descriptions reference 1D signals, the results can be extended to higher dimensions (i.e. 2D image filters).


what are the applications of fourier transform in modulation what are the applications of Fourier transform in modulation ? What are modulation schemes that use Fourier transform? I guess, FSK & PSK are related to Fourier transform since it both have frequency and phase in their names respectively.

OFDM is usually modulated by using an IFFT (or FFT), using the DFT basis vectors for the multi-carrier orthogonality.


Question on covariance matrix of 2 spatial signals Every time I think I have understood the covariance matrix, someone else comes up wih a different formulation.

I am currently reading this paper:

  J. Benesty, "Adaptive eigenvalue decomposition algorithm for passive acoustic source localization", J. Acoust. Soc. Am. Volume 107, Issue 1, pp. 384-391 (2000)

and I have come across a formulation I do not quite understand. Here, the author is constructing the covariance matrix between two signals, $x_1$, and $x_2$. Those two signals are from different sensors.

For the covariance matrix of one signal, I know that we can get it by calculating the regression matrix, and then multiply it by the Hermitian of that same matrix, and dividing by $N$, the length of the original vector. The size of the covariance matrix here can be arbitrary, with maximum size being $N\times N$.

For the covariance matrix of two spatial signals, if we place the first signal in the first row, and the second signal in the second row of a matrix, then multiply by its Hermitian, and also divide by $N$, then we get a $2\times 2$ covariance matrix of both spatial signals.

However, in this paper, the author computes what looks like four matricies, $R_1{}_1, R_1{}_2, R_2{}_1$, and $R_2{}_2$, and then puts them into a super matrix and calls that the covariance matrix.

Why is this so? Here is an image of the text:

If you have two signal vectors $x_1[n]$ and $x_2[n]$ each of $N$ elements, then there are two different things we can consider.

How do the quantities $\displaystyle\sum_{n=1}^N x_i[n]x_j[n],~ i, j \in \{1,2\}$ compare?  In particular, when the signals are noisy and the noises can be considered
to be jointly stationary (or jointly wide-sense stationary), these quantities can be used to estimate the noise variances in the two signals as well as the covariance of the noises at any fixed sampling time.  This is what you get from the $2\times 2$
covariance matrix
$$R_{2\times 2} = \left[\begin{matrix} \sigma_1^2 & C\\
C & \sigma_2^2\end{matrix}\right].$$
The noise in $x_1[n]$ has variance
$\sigma_1^2 = R_{1,1}$ which might be different
from $R_{2,2} = \sigma_2^2$, the variance of the noise in $x_2[n]$.  But the noises are
correlated with covariance $R_{1.2}=R_{2,1} = C$.  Now, if we plan on doing things
with just what happens at $n$, ignoring whatever might be happening at $n-1$ or $n+1$
etc., then this is all the information we need.
Unless the noise is known to be (or assumed to be) white noise so that noise
samples from different sampling instants are independent (and therefore uncorrelated)
or we simply assume uncorrelated noise samples, there is information that we
are ignoring by not considering the correlation between $x_1[n]$ and $x_1[m]$,
samples from the same process at different times or locations, and the correlation between $x_1[n]$ and $x_2[m]$, samples from the two processes at different times
or locations.  This additional information might lead to a better estimate/solution.
Now  we have a total of $2N$ noise samples,
and therefore a $2N \times 2N$ covariance matrix to consider.  If we
arrange matters the way the authors did, we have
$R_{\text{full}} = E[XX^T]$ where
$$X = (x_1[1],x_1[2],\ldots,x_1[N],x_2[1],x_2[2],\ldots,x_2[N])^T
= (\mathbf x_1,\mathbf x_2)^T$$
and so
$$R_{\text{full}} = \left[
\begin{matrix}R_{\mathbf x_1,\mathbf x_1} & R_{\mathbf x_1,\mathbf x_2}\\
R_{\mathbf x_2,\mathbf x_1} & R_{\mathbf x_2,\mathbf x_2}\end{matrix}
\right]$$ where $R_{\mathbf x_i,\mathbf x_j} = E[\mathbf x_i\mathbf x_j^T]$.
Note that $R_{\mathbf x_i,\mathbf x_j}$ is, in essence, the
cross-correlation function of $(x_i[1],x_i[2],\ldots,x_i[N])$
and $(x_j[1],x_j[2],\ldots,x_j[N])$ if $i \neq j$ and the
autocorrelation function if $i = j$.  If the noise processes are white and
uncorrelated except when $n = m$, then
$$R_{\text{full}} \rightarrow R_{\text{simple}} = \left[ \begin{matrix}\sigma_1^2I & CI\\
CI & \sigma_2^2I \end{matrix}\right]$$
where $I$ is the $N\times N$ identity matrix, and
$\sigma_1^2, \sigma_2^2$ and $C$ are as defined in Item 1 above.
How realistic this noise model might be is something for the end user
to determine. If the model is realistic, then nothing is gained
by looking at the $2N\times 2N$ matrix $R_{\text{full}}$
since all the information is there
in the $2\times 2$ matrix $R_{2\times 2}$ of Item 1 above.  Ditto if the model is
unrealistic but we don't intend to (or are unable to)
use all the information in the full
$2N\times 2N$ matrix $R_{\text{full}}$; we will make do
with just $\sigma_1^2, \sigma_2^2$ and $C$ of Part 1
for which we don't need $R_{\text{full}}$ or $R_{\text{simple}}$,
just $R_{2\times 2}$.


What is the beneficient of taking the FFT of an ecg signal? Our teacher gave us a homework.
He says first find an ecg signal, then take FFT of that signal?

Why do we need to take FFT of an ecg signal? Is it that we can cancel noise only if we know the frequency spectrum of the signal?

I'm not a doctor, I'm an engineer. However, I think I see the point of this.

An electrocardiograph is used as a diagnostic tool to assess the state of the heart and cardiovascular system. Several signals are measured, each of which can be used to make judgements about the health of the patient. Most of the signals measured in an ECG are periodic (in other words, they occur at a given frequency). The DFT/FFT is basically a way of analyzing the frequency content of a signal. So, it follows that by taking the DFT/FFT of an ECG, you can easily see those frequency components.

From the Wikipedia article on ECGs, here are some common components of an ECG and their associated frequencies:

P Wave: 0.08s, 12.5Hz

QRS Complex 0.09s, 11Hz

T Wave: 0.16s, 6.25Hz

I don't actually know the medical meaning of those features, but in any case, they're periodic, meaning you can analyze them using spectral analysis.


Why is the signum function $=2u(t) - 1$? Basically I couldnt find proof for this anywhere though it's a very simple and basic equation.

The signum function is defined by

$$\text{sgn}(t)=\begin{cases}-1,&t<0\\0,&t=0\\1,&t>0\end{cases}$$

Using the half-maximum convention, the unit step function is defined by

$$u(t)=\begin{cases}0,&t<0\\\frac12,&t=0\\1,&t>0\end{cases}$$

From these two definitions it should be obvious that

$$\text{sgn}(t)=2u(t)-1$$

must hold.


Using all-pole filter to model the Room Impulse Response Why is all pole model pretty useful in modelling room acoustics?

Is it related to reverberation?

It's not because of reverberation.

When you want to model the Frequency Response of the room, it's common to simplify your approximation by using either all-pole or all-zero models. You don't want to use the full zero-pole model.

To get some intuition:

zeros correspond to time delays and antiresonances
poles correspond to resonances of your Room Response

In practice all-zero models are not being used due to various reasons, such as:

required filter length is comparable to the IR length, and almost 40x the length of corresponding all-pole filter
filter will be valid only for specific distances and positions between the source and receiver (remember: time-delays).

That is why the all-pole model is used instead. As mentioned above, poles correspond to the resonances, i.e. standing waves, which are:

independent of the source location (quite intuitive)
independent of the receiver location (except of the nodes)

Additionally the required filter length is way less than in case of all-zero models. According to Mourjopoulos, for $RT_{60}\approx0.5 \mathrm{s}$, the required order is within the range of $50 < N < 500$.  The same author, concludes that all-pole filters are easier to manipulate than all-zero filters, due to their filter length. Morevover, author mentioned that all-pole filters are sufficient approximation than using raw impulse response data.

Here is some literature:

  Mourjopoulos J. - On the Variation and Invertibility of Room Impulse
  Response
  Functions

  Mourjopoulos J., Paraskevas M. A. - Pole and Zero Modelling of Room
  Transfer
  Functions


How does companding compress dynamic range? Regarding when you 'compress' the signal, I typically see diagrams like this:

And the companding curve:

I can only see how the companding curve achieves a higher average SQNR by rounding down higher amplitudes more than lower amplitudes when making the quantization, according to a normal PDF, making the quantisation more efficient. I don't see how it reduces the dynamic range of the signal to the 'compressed dynamic range'. The only dynamic range compression I can see is that when the higher amplitudes are expanded again, they'll be expanded to the figure they were rounded down to during quantisation, reducing the upper dynamic range, but not the lower.
The dynamic range compression appears to be referring to the dynamic range of the signal after the compressor and before the expander, but this doesn't make sense. How can PCM samples have a dynamic range that isn't anything you want it to be.. and as the samples are uniform, surely a high rounded down amplitude now has a higher relative amplitude (sample 255/256 for instance whereas it was originally rounded down much more), and of course a lower amplitude is now a much higher amplitude.
What confused me is that there is a correspondence between 1.0 input amplitude and 1.0 output amplitude, which doesn't compress the dynamic range at all.
I hope that the key word in the output is 'relative' i.e. relative to the output range and it's actually mapping it onto a smaller dynamic range, though it looks like it's the same range. It would make sense to me if the dynamic range compression comes from the dynamic range that the samples are being mapped to is smaller than the input, and the companding curve is just a way to do this more efficiently than linear.
Why is the companding curve even called 'compression', surely the compression part is mapping it to a smaller dynamic range, or when using digital compression, the compression part is the reduction of bits in the PCM sample, and you can compress with linear companding, which wouldn't be considered compression at all if you were defining compression to be non uniform quantisation.

Companding is used when you want a signal to noise ratio (SNR) that's largely independent of the signal level.
Quantization noise is constant, so if the signal level varies, so does the SNR.
Let's consider an audio signal that's quantized at 8 bit. The quantization noise sits at -53dBFS. If you have a loud signal, i.e. -10dBFS you have a healthy 43 dB SNR. However, for a soft signal at -40dBFS, the SNR is very poor at 13 dB.
If you compand you can get an SNR of maybe 30dB for both cases, which is much better. The audibility of the quantization noise is NOT a function of the absolute noise level but the noise level relative to the signal. You can tolerate more noise, if the signal itself is louder since the louder signal masks more of the noise.


Offset in Carrier and Timing Is my understanding about carrier and timing offset right?

Suppose my transmit  frequency is 100 MHz and if my reference is 1ppm; then the transmit frequency has +/-100Hz offset.

Similarly, if my receive frequency is 100 MHz and if my reference is 2 ppm; the receive frequency  has +/-200Hz offset.

So the carrier frequency offset is +/- 300 Hz, i.e 600 Hz; does this have to be recovered in receiver carrier frequency recovery?

Similarly, if the symbol rate is 100KHz with reference 1 ppm then the transmitter has +/- 0.1 Hz symbol offset and the receive symbol rate is 100kHz with reference 2 ppm it has +/- 0.2 Hz symbol offset.

So the symbol frequency offset is +/-0.3 Hz, i.e 0.6 Hz. Is this symbol recovery? What is symbol phase recovery? Is there any effect of sampling frequency offset also?

The calculations the OP is using is correct to predict the maximum carrier and symbol timing offset due to reference clock offsets in the transmitter and receiver.

For sampling clock offsets (both in the DAC in the transmitter and the ADC in the receiver), any offset of the sampling clock also introduces a carrier frequency offset that will be in proportion to the carrier that the analog signal is at when sampled according to $f_a/f_s$, where $f_a$ is the analog carrier frequency (of offset from $f=0$) and $f_s$ is the sampling rate. Thus a ppm change in the sampling clock will result in the same ppm change of the carrier of the sampled signal, given the actual carrier the signal is at when sampled. For example if an analog signal was translated down to a 25 MHz IF (meaning a 25 MHz carrier), and was sampled with 100 MHz sampling clock, then $f_s=100e6$ and $f_a=25e6$ in this case, and the fractional frequency of the digital IF would be $f_a/f_s = 0.25$ samples/cycle.  If the sampling clock was actually 1% higher, or 101 MHz, the fractional frequency would 1% lower given $f_a/f_s \approx 0.2475$. Be very careful to reference the actual analog carrier and not the digital IF which can be different in undersampling applications. For example, if the signal was at a 125 MHz analog IF and sampled with a 100 MHz clock to produce a 25 MHz digital IF, the fractional frequency with no clock error would be $(f_a-f_s)/f_s = 0.25$ but if the sampling clock was 1% higher, or 101 MHz, the fractional frequency would be $(f_a-f_s)/f_s \approx 0.2376$ which is 5% lower! This works out to to be the clock error in ppm multiplied by the ratio of the analog IF frequency to it's equivalent frequency in the first Nyquist zone, in this case $125/25 = 5$.

Additional carrier frequency offset can be introduced through channel effects such as Doppler due to the motion of the transmitter and receiver that also must be taken into consideration for determining the maximum carrier and symbol timing offsets.

Once at baseband, a symbol phase offset would typically be referring to the carrier phase offset itself (rotation of the symbol) since an added phase at the carrier directly results in an added phase of the symbol. Even though it would be a referred to as a carrier phase offset, it may not have originated at the original carrier of the signal but can be introduced by other factors in the processing of the signal. If this carrier phase is changing with time, then this is a carrier frequency offset given frequency is the derivative of phase with time. Symbol timing offset at baseband in contrast would refer to the offset between the ideal sampling location for minimum error rate versus the actual sampling location. (literally a sampling time offset, or delay), and there is no rotation involved when the sampling location on the baseband signal is delayed in time. Similarly if a waveform was sampled directly at baseband, only timing offset would be introduced if there were phase or frequency offsets in the sampling clock since the IF frequency at baseband = 0 (and a frequency offset in the sampling clock would generate a varying time offset on the baseband constellation but not a rotation). Any small residual carrier offsets would be affected just like the IF frequencies above if they were then sampled with a clock that has offset errors, but this will typically be negligible given the effective IF frequency is so small.


The Effect of the Standard Deviation ($ \sigma $) of a Gaussian Kernel when Smoothing a Gradients Image I am trying to understand why smoothing an image with a gaussian kernel of different sigma values and then computing the gradients of the smoothed image leads to "thicker" trails. In the image below, the left image is produced by convolving the image with the derivative of a gaussian kernel where $\sigma = 1$ while the right picture shows image gradients when the image is convolved with a gaussian kernel where a $\sigma = 2$.
From what I understand, the gaussian kernel with a higher $\sigma$ value tends to make pixels look "more" like it's neighbors, because it places higher weights on pixels that are further away from the middle.

Let's analyze it in 1D as the intuition is the same.
First, let's have a look on a few different Gaussian Kernels:

As expected, they are wider as the Standard Deviation (STD) increase.
It means that when the kernel is applied using the convolution, more information is aggregates from farther samples. On the other side it means data is spread.
Now, in your images a gradient is a bump. So we can approximate it by a step function:

Let's have a look at the output of the convolution of the kernels from above with the step function:

As we can see, the higher the STD the data is more spread, Though the amount of energy is the same.
In image we expect it to be seen as data is more blurred, wider, while being darker (As the values are lower).
The full code is available on my StackExchange Signal Processing Q70725 GitHub Repository.


Obtaining databases for speech emotion recognition analysis For a research work on speech emotion recognition, I need the database of CASIA and EESDB for simulations. I could not find those; how can I get them? Thanks.

CASIA can be obtained from Chinese LDC:

http://shachi.org/resources/27

Overall, it is easier to get experiments going with EMO-DB or to collect your own data from Youtube. Unsupervised learning from videos, like youtube videos tagged by uses with emotional tags might make more sense in a long term.


Apply Low pass Butterworth filter in Python I want to use a low pass Butterworth filter on my data but on applying the filter I don't get the intended signal. Here is the dummy code:

Signal A:

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
a = np.linspace(0,1,1000)
signala = np.sin(2*np.pi*100*a) # with frequency of 100
plt.plot(signala)

Signal B:

signalb = np.sin(2*np.pi*20*a) # frequency 20
plt.plot(signalb)

Let's combine signal A and B now to get signal C

signalc = signala + signalb
plt.plot(signalc)

The resultant signal C looks like this

Let's now apply the filter:

b, a = signal.butter(5, 30, 'low', analog = True) #first parameter is signal order and the second one refers to frequenc limit. I set limit 30 so that I can see only below 30 frequency signal component
output = signal.filtfilt(b, a, signalc)
plt.plot(output)

On applying above butter filter, I get an empty plot as

I don't understand where I am doing wrong? Any  help will be appreciated.

You should not be using the analog filter - use a digital filter instead. You want the filter to be defined in Z-domain, not S-domain. Also, you should define the time vector with known sampling frequency to avoid any confusion.

The design of the digital filter requires cut-off frequency to be normalized by fs/2.

Here is a working example:

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

fs = 1000  # Sampling frequency
# Generate the time vector properly
t = np.arange(1000) / fs
signala = np.sin(2*np.pi*100*t) # with frequency of 100
plt.plot(t, signala, label='a')

signalb = np.sin(2*np.pi*20*t) # frequency 20
plt.plot(t, signalb, label='b')

signalc = signala + signalb
plt.plot(t, signalc, label='c')

fc = 30  # Cut-off frequency of the filter
w = fc / (fs / 2) # Normalize the frequency
b, a = signal.butter(5, w, 'low')
output = signal.filtfilt(b, a, signalc)
plt.plot(t, output, label='filtered')
plt.legend()
plt.show()


is $F(x,n)=x(-n + 2)$ time invariant? My book with answers to proakis and manolakis book says that $F(x,n)=x(-n + 2)$ is time invariant while my calculation shows that it is time variant. I have repeated them many times and I am mostly convinced that book is wrong. Can somebody please check my calculations?

First we shift output signal:

$$F(x, n - k) = x(-(n-k) + 2) = x(-n+k+2)$$

now lets put $x'(n) = x(n-k)$ and check what happen when we input shifted signal:

$$F(x', n) = x'(-n+2)=x(-n+2-k)=x(-n-k+2)$$

since $F(x, n - k) \neq F(x', n)$ this system is time variant. Is my method correct?

Your reasoning is correct and we arrive at the same conclusion.

The system's output is given by

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

Let

$$x_2[n]=x[n-k]\tag{2}$$

The corresponding output is

$$y_2[n]=x_2[-n+2]=x[-n-k+2]\tag{3}$$

From $(1)$ the shifted output is

$$y[n-k]=x[-(n-k)+2]=x[-n+k+2]\tag{4}$$

Since $(3)$ and $(4)$ are not identical, the system is time-varying.


Digital filters with more zeros than poles I am having trouble wrapping my head around digital filters with different orders of numerator and denominator.  Let me know if any of these points is wrong:

All (digital or analog) transfer functions have the same number of poles and zeros, if you include the ones at infinity.  So $H(s) = 1/s$ has a pole at the origin and a zero at infinity, which is important for visualizing the amplitude response surface in the S plane.
But usually when we say "number of poles" or zeros we mean finite poles or zeros.  So $H(s) = 1/s$ is considered to have one pole and no zeros.
Digital filters designed by bilinear transform from analog filters always have the same number of poles and zeros (and none are ever at infinity).  But digital filters can also be made with different orders of numerator and denominator.
To find the poles and zeros of a digital transfer function in the general case, you first must express it as positive powers of z ("controls engineer format") and then find the roots of the numerator and denominator, same way as an analog filter.  So for example, a single-sample delay in "DSP engineer format" $H(z) = z^{−1}$ is rewritten as $H(z) = 1/z$, which shows that it has a pole at the origin and a zero at infinity, so it's considered to have "one pole and no (finite) zeros".

So then I become confused:

FIR filters are described as "all-zero filters". 1 2 They can be represented as a transfer function like $$H(z) = b_0 + b_1 z^{-1} + b_2 z^{-2}$$  But if you convert to positive powers of z to find the poles and zeros, you get: $$H(z) = \frac{b_0 {z}^{2} + b_1 z + b_2}{z^2}$$ which has just as many poles at the origin as there are zeros.  What is the significance of these poles?  These are each the $H(z) = z^{−1}$ delay elements used to produce the feedforward signals?  Seems like the poles end up at the origin if the delay elements are not fed back to the input?  So FIR filters are not actually all-zero filters?
Similarly, maxflat filters (Selesnick-Burrus generalized Butterworth) are described as having "more zeros than poles", which is supposed to be computationally advantageous.   (Why?)  The Matlab example produces b = [0.0950    0.2849    0.2849    0.0950] and a = [1.0000   -0.2402].  I think this is "negative powers of z" format, so this would represent a transfer function $$H(z) = \frac
{0.0950 + 0.2849 z^{-1} + 0.2849 z^{-2} + 0.0950 z^{-3}}
{1.0000 - 0.2402 z^{-1}}
$$
But again, if you convert to positive powers of z, you get: $$H(z) = \frac
{0.0950 z^3 + 0.2849 z^2 + 0.2849 z + 0.0950}
{1.0000 z^3 - 0.2402 z^2}
$$ which has 3 poles again, 2 of which are at the origin.  So the number of poles and zeros is again the same, and they're all finite, too.

If adding a zero to a digital transfer function always also adds a finite pole, I don't understand how a filter with "more zeros than poles" can exist, or be advantageous.  Might as well use those poles to affect the frequency response if you have them?

Is there some convention where poles at the origin are not included in the tally, like the way poles at infinity are not included?  If so, why?

If you consider the transfer function of a causal IIR filter

$$H(z)=\frac{B(z)}{A(z)}=\frac{\sum_{m=0}^M b_mz^{-m}}{\sum_{n=0}^N a_nz^{-n}},\quad a_0=1$$

then you always get the same number of poles and zeros, regardless of the choice of $M$ and $N$ (as already pointed out by Robert). However, what is meant by a system with "more zeros than poles", is a system with a numerator degree $M$ greater than the degree $N$ of the denominator. In this case "poles" refers to the poles away from the origin. It is only these poles that have to be implemented. So a system with some poles at the origin is cheaper to implement than a system with all its poles away from the origin. So your argument that you "might as well use those poles ... if you have them" does not really hold, because poles at the origin come without any implementation costs.

The remaining question is why such filters with more zeros than poles (away from the origin) can be useful. Let me just give two examples:

In implementation as well as in the design process, poles can give a lot of trouble (instability, noise enhancement, local optima in the design process, etc.). On the other hand, FIR filters with good frequency selective properties must often have a very large degree (resulting in a large delay and in high implementation costs). Adding just a few poles can give us the best of both worlds: better filter behavior at a greatly reduced (numerator) degree $M$, and only little trouble because there are only very few poles.
For frequency selective filters with a small desired phase distortion in the passbands, IIR systems with $M>N$ are very useful, because the poles ideally contribute only to the passbands, whereas the zeros must contribute to the passbands (for phase equalization) as well as to the stopbands (zeros on or close to the unit circle). Consequently, we need poles only at angles corresponding to passband frequencies, but we need zeros everywhere.

The figure below shows such a system (lowpass filter with approximately linear passband phase response, $M=12, N=6$):


Determining if the system is linear and time invariant The Fourier transforms of input $x[n]$ and output $y[n]$ of a discrete-time system are related by the following equation: $$Y(e^{j\omega}) + (e^{-j\omega})Y(e^{j\omega})= X(e^{j\omega}) - X^{'}(e^{j\omega})$$ Determine if the system is (i) linear (ii) time-invariant.

Hint : Use following properties of Discrete-Time Fourier Transform :

If $DTFT\{x[n]\} = X(e^{j\omega})$, then $DTFT\{x[n - n_0]\} = e^{-j\omega n_0}X(e^{j\omega})$
$IDTFT\{j\frac{\partial{}}{\partial\omega}X(e^{j\omega})\} = n x[n]$


Calculating frequency shift of short bursts I am attempting to figure out the best method of determining the frequency shift of small burst signals. I have an incoming signal which is 1Mhz, and only in 0.0001s windows, seeing 10kHz shifts, which I am trying to measure using an FPGA and ADC. I have tried to run the signal through a BPF and FFT but the computing time is long enough that I am missing a bunch of signals and losing accuracy. Now I'm trying to take the raw signal and measure phase change, and/or frequency by count the crossings of a threshold value. I've tried looking at larger windows, but the bursts are small enough that they essentially get erased.

Are there high-speed methods for measuring signals like this? I'm afraid too much decimating is required for the FFT

Are there ways to measure phase/frequency that are a bit more flexible? My resolution is now 10kHz and only ~0.3% of samples have shift. Better resolution is erasing the shift since windows are longer.

TIA for any comments!

Try using the instantaneous frequency

$$f [n] =\dfrac {F_s/2} {\pi}Arg (x [n] \cdot x^*[n-1]) $$

and average it over the burst.


Question on understanding the time-bandwidth of a signal I read from a paper (Ozaktas, H.M., et al. (1996) Digital Computation of the Fractional Fourier Transform) that says
The time-bandwidth product can be crudely defined as the product of the temporal extent of the signal and its (double-sided) bandwidth. It is equal to the number of degrees of freedom and the number of
complex numbers required to uniquely characterize the signal among others of the same time-bandwidth product.
I understand the first statement, but the second statement is hard to interpretate.
First, why is the time-bandwidth product equal to the number of degrees of freedom and what is the number of degrees of freedom of a signal? Does it hold for both the deterministic and random signal?
Second, why is there a (minimum) number of complex numbers required to characterize a signal? What if the number is not reached?
Hope you can help clarify these confusions.

This is not a rigorous answer.

First, why is the time-bandwidth product equal to the number of
degrees of freedom and what is the number of degrees of freedom of a
signal? Does it hold for both the deterministic and random signal?

The degree of freedom is the number of samples needed to represent a complex continuous-time signal. The Nyquist Shannon sampling theorem states that for any signal confined within $\left[-\frac{\Delta f}{2},+\frac{\Delta f}{2}\right]$, we need at least $\Delta f$ sample per time unit to perfectly reconstruct it.
If in the time domain, the signal can be approximated by its duration $\left[-\frac{\Delta t}{2},+\frac{\Delta t}{2}\right]$, the number of samples needed is thus $\Delta f \Delta t$, which is the time-bandwidth product defined in [Ozaktas et al.].

Second, why is there a (minimum) number of complex numbers required to
characterize a signal? What if the number is not reached?

If the number is not reached, the discrete-time signal represented by the insufficient number of samples cannot be used to reconstruct the continuous-time signal. An aliasing will be observed, for example.


Amplitude of a signal. [-1, 1] vs variants and MFCC implications I've seen two sort of audio signals, one where the y-axis takes values between -1 and +1. I'd think this means 1 means loudest? Not really sure

And the second one I saw goes between $-10 000$ and $+10 000$. What do these values actually mean? Are they arbitrary?

I'm doing work on audio features using MFCC but it seems that when I take a file with values between $[-10k,+10k]$, I get normal cepstral coeffecients of which we take the first $12$. But when I take audio with values in $[-1, 1]$, the cepstral coeffecients are in a negative range.

What do these values actually mean? Are they arbitrary?

Yes, pretty much: they will just be the scale of numbers your system works with. In floating point systems, we often see that samples get normalized to [-1,+1], whereas in fixed-point system, it's often things like [-2⁻¹⁵,+2¹⁵-1], depending on the bit width of the samples to begin with.

So, this is purely conventional. Most algorithms don't care about the absolute scale of things.

So, if your algorithm does depend on amplitudes, well, normalize your signal before you process it.


Metal Defect Detection I am developing a machine vision project to find defects on metal surfaces. I want to learn and implement algorithms which commercial systems are use for metal surface defect. My sample object as below:

I am using OpenCV. I could find sharp defects with standard edge detection methods. But it is not sufficient. There is some samples of commercial solutions:
- http://www.keyence.co.uk/products/vision/vision-sys/cv-5000/features/movie06.jsp
- http://www.keyence.com/products/vision/vision-sys/cv-5000/features/movie09.jsp

Or like this picture:

Do you have any idea about algorithms which is used in commercial systems?

As the video says, they use the classical lowpass filtering for texture smoothing, then shading correction to lessen the effect of uneven illumination, followed by binarization for defect segmentation.

A completely different technique is used for the bearing: edge detection and (presumably) assessment of deviations from a smooth line.

The scratch detection on the heavily textured surface is much more challenging and I am not sure that a commercial system can handle it.

For the best of my knowledge, the technique in "Fast Detection of Curved Edges at Low SNR, Nati Ofir, Meirav Galun, Boaz Nadler and Ronen Basri" could do.


BER performance of fractional equalizers for colored noise I am trying to compare the BER performance of fractional and symbol spaced equalizers for channel equalization in  colored noise scenario. I have found from simulation that for positively correlated noise, fractional equalizers perform better than symbol spaced equalizers. It is probably due to fact that with more samples we have more information in case of fractional equalizers. Is there any way of theoretically prove that fractional equalizers are better than symbol spaced equalizers for colored noise?
Thanks

Yes of course,  that's true.

FS equalizer sample the signal in rate faster than Nyquist rate which provide more information about the signal. It's similar when you  increase the sampling frequency, you get better results.

If you need to understand that idea, you can read about oversampling, fractional sampling, and if you are working on sparse channel read you compressive sending working method.


Correlation between frequency and initial phase estimates using FFT Why is there a high negative correlation between frequency estimation errors and initial phase estimation errors when measured using FFT?
I have a simple code showing this negative correlation.
is there a known explanation for that?
here is the code

clear all;clc;close all;
%%
T = 0.2;
fs = 500;
f = 200;
ts = 1/fs;
snr = 5;
n = fs*T;
t = linspace(0,T-ts,n);
theta = 0.1;
y = exp(1i*2*pi*(f.*t + theta));
for count = 1:2000
    count;
    yy  =  awgn(y,snr,'measured');
    NFFT = 2^16;
    fVals = fs/NFFT.*[0:NFFT-1];
    X = fft(yy,NFFT)./n;
    plot(abs(X))
    [~,M] = findpeaks(abs(X),'MinPeakheight',0.8);
    measuredAngle = angle(X(M))/2/pi;
    errfreq(count) = fVals(M) - f;
    errph(count) = measuredAngle - theta;
end
cov(errfreq,errph)
corrcoef(errfreq,errph)

ans =
    1.0000   -0.8506
   -0.8506    1.0000

It's a spurious artifact due to not referencing your phase generation to the center of your data window, and not doing an fftshift before the fft to place your phase reference point at fft input index 0.

e.g. For varying frequencies, either you have to know, generate, or want the phase at location T/2 of your continuous sinusoid.  The phase result of an FFT after an fftshift will produce an estimate of the phase at T/2 for comparison.  For continuous signals, when using this fftShift, interpolation to improve phase estimation referenced to the center point is possible (Sinc kernel better than linear or polynomial).

This has to be done anytime the input signal isn't absolutely locked to be perfectly integer periodic within the fft length in order to get reasonable phase measurement results.  e.g. for any sweeping, variable, or unknown frequencies.

Otherwise the phase will toggle during the frequency sweep.  This is because non-integer periodic sinusoids have a circular discontinuity at the window boundaries, and this discontinuity changes sign during frequency sweeps, causing the FFT phase result to have a discontinuity.  An fftshift moves this discontinuity away from the beginning or phase measurement point of an FFT, which is the reference point where both the cosine and sine, or real and imaginary components of the FFT basis vectors start with a phase of zero.


Difference Fourier Shift Theorem and Cross Correlation Theorem can someone explain me the difference between the Fourier Shift Theorem and the Cross Correlation Theorem?

  The shift theorem says that a delay in the time domain corresponds to a linear phase term in the frequency domain

source

So in case of image processing I can find out the shift of two images, which are shifted against each other. I only need to calculate the normalized cross power spectrum of this two images and do an invers Fourier Transform to get the peak, which represents the shift of these to images.

The Cross correlation Theorem says

  the cross-correlation between two signals is equal to the product of fourier transform of one signal multiplied by complex conjugate of fourier transform of another signal.

source

So what is the difference between these two theorems?

With both I can for example calculate the shift between two images. I can do it, in both ways, with the Phase Correlation (cross-power-spectrum), right?

First of all, what you quote as the Cross correlation Theorem

  the cross-correlation between two signals is equal to the product of fourier transform of one signal multiplied by complex conjugate of fourier transform of another signal

doesn't make much sense: the product mentioned above gives the _cross_power spectral density function, and not the cross-correlation function at all. The cross-correlation function is the inverse Fourier transform of the cross-power spectral density function, namely, the product mentioned above.

That being said, the key difference between what you are doing and what the Cross-Correlation Theorem is that you are using a special case that is perhaps better termed the Autocorrelation function theorem.  Notice the use of the word another (next to last word in the highlighted wording above). When another is replaced by the same, we get that the autocorrelation function is the inverse Fourier transform of the power spectral density.


Convert a FIR to an equivalent IIR Is there a way to convert a FIR to an IIR filter with the most similar behavior?

I would say that the answer to your question - if taken literally - is 'no', there is no general way to simply convert an FIR filter to an IIR filter.

I agree with RBJ that one way to approach the problem is to look at the FIR filter's impulse response and use a time domain method (such as Prony's method) to approximate that impulse response by an IIR filter.

If you start from the frequency response then you have lots of methods for designing IIR filters. Even though it was published about 25 years ago, I believe that the method by Chen and Parks is still one of the better ways to approach the design problem. Another very simple method for the frequency domain design of IIR filters is the equation error method, which is described in the book Digital Filter Design by Parks and Burrus. I've explained it in this answer.

If the phase response is of importance to you, then one problem you will be facing when designing IIR filters in the frequency domain is the exact choice of the desired phase response. If the overall shape of the desired phase is given you still have one degree of freedom, which is the delay. E.g., if the desired phase is $\phi_D(\omega)$, and the desired magnitude is $M_D(\omega)$ then your desired frequency response can be chosen as

$$H_D(\omega)=M_D(\omega)e^{j(\phi(\omega)-\omega\tau)}\tag{1}$$

where $\tau$ is an unknown delay parameter. Of course you can say that if $\phi_D(\omega)$ is given then you don't want to modify it with an additional (positive or negative) delay. But it turns out that in practice the average delay is not always important, and - more importantly - for certain values of $\tau$ your approximation will be much better for a given filter order than for others. So the delay $\tau$ can become an additional design parameter and should be chosen optimally or at least reasonably.

I've written a thesis on the design of digital filters with prescribed magnitude and phase responses. One chapter deals with the frequency domain design of IIR filters. That method can be used to design IIR filters with approximately linear phase in the pass-bands, or to approximate any other desired phase (and magnitude) response. The filters are not only guaranteed to be stable, but you can also prescribe a maximum pole radius, i.e., you can define a certain stability margin. You can also find this method in a paper published in the IEEE Transactions on Signal Processing.


Goertzel algorithm magnitude I'm reading a Goertzel algorithm paper and an embedded article. I have several questions:

I don't understand why the the magnitude is calculated as follows:

$$\text{Magnitude} = y^2[n-2] + y^2[n-1] - a * y[n-2]y[n-1]$$
where $a$ is some coefficient. Can someone elaborate on why this is the expression?

Minimum block size If I have a sampling rate of 4000 Hz and target frequency of 50 Hz, what is the minimum block size?
If block size more that magnitude will more accurate right?

[1] It's a notch filter. It attenuates everything but the frequency component you're interested in. The output is summed up over some time to make the tone detection more reliable.
[2] There's no minimum block size. There's a requirement about the amount of time spent summing up the output. What you'll see is that the summed output over say 205 samples (at fs 8000kHz), grows much more if the input contains the frequency component you're interested in detecting.