Downsampling: Mathematical derivation The problem I am having is related to sample rate conversion and more precise to sample rate reduction. I have been working on the paper Interpolation and Decimation of Digital Signals Tutorial Review in [1] and A digital signal processing approach to interpolation in [2] where the maths make my understanding difficult.

To start with,bellow is defined a new sequence $w'(n)$

$$
w'(n)=
\begin{cases}
w(n), &\text{if}\quad n=0,\pm M, \pm 2M, \ldots \\[2ex]
0, & \text{otherwise}
\end{cases}\tag{1}
$$

that is non zero only at multiples of $M$ of the old sampled sequence. Then it is stated that a convinient representation of $w'(n)$ is:

$$
w'(n)=w(n)\left(\frac{1}{M}\sum^{M-1}_{l\ =\ 0}e^\frac{j2\pi ln}{M}\right)\tag{2}
$$

The term in the brackets corresponds to a discrete Fourier series representation of a periodic impulse train with a period of $M$ samples. That in my understanding means that is a multiplication by a $\delta(n)$ in time domain if the term in the brackets refers to the inverse DFT of one $(1)$. If that is the case then shouldn't that be a convolution?

Then an expression for elaborating on the $z$-domain is derived:

$$
 y(m)=w'(Mm)=w(Mm)
$$

Then $\mathcal Z$-transform is utilized after:

\begin{align}
Y(z)&=\sum^{+\infty}_{m=-\infty}y(m)z^{-m}\\
&=\sum^{+\infty}_{m=-\infty}w'(Mm)z^{-m}\\
&=\sum^{+\infty}_{m=-\infty}w(m)\left(\frac{1}{M}\sum^{M-1}_{l\ =\ 0}e^\frac{j2\pi lm}{M}\right)z^{-\frac{m}{M}}\\
&=\frac{1}{M}\sum^{M-1}_{l\ =\ 0}\left(\sum^{+\infty}_{m=-\infty}e^\frac{j2\pi lm}{M}z^{-\frac{m}{M}}\right)
\end{align}

and we end up here.

How to conclude the result bellow?
What are the steps missing?
Where has $m$ gone?
And why do we have the exponential as argument of $W$ in equations $(3)$?

$$
Y(z)=\frac{1}{M}  \sum^{M-1}_{l\ =\ 0}  W\left(e^\frac{j2\pi l}{M}z^{\frac{1}{M}}\right)\tag{3}
$$

[1] R.E.Crochiere and L.R.Rabiner,"Interpolation  and Decimation of Digital Signals Tutorial Review",Proceedings of the IEEE,vol.69,pp 300-331,March 1981

[2] R.W.Schafer and L.R.Rabiner,"A digital signal processing approach to interpolation",Proceedings of the IEEE,vol 61,pp 692-702,June 1973

@MattL 's link solves your problem, but it seems you need more help. Then I would like to provide here the same solution in different approach, expecting it might help.
Analysis of a downsampler block (more specifically the decimator or the compressor block) can be done in a number of ways. I prefer the following approach which I find it easier.
The first step of the analysis begins by observing that the structure of the figure a and figure b produce exactly the same outputs $y[n]$  given the same inputs $x[n]$, hence they are equivalent:
$$ \boxed{ x[n] \rightarrow \boxed{\downarrow M} \rightarrow y[n] ~~~~~~,~~~~~~~~~~~ \text{figure-a} }$$
$$ \boxed{
x[n] \rightarrow \boxed{ (\times)} \rightarrow v[n] \rightarrow \boxed{\downarrow M} \rightarrow y[n] ~,~ \text{figure-b} \\
  ~~~~~~~~~~~~~~~\uparrow w[n] }$$
Where the multiplier $w[n]$ is a periodic impulse train with a period of $M$ samples, which constitudes a discrete-time counterpart of a continuous-time impulse train sampler.
$$ w[n] = \sum_k \delta [n- kM] $$
Such kind of a structural modification, enables an unambiguous mathematical analysis at the expense of an increased number of steps and structural complexity. We shall analyse the figure-b and map its result to figure-a which is what we want indeed.
Then, as the second step, express the periodic impulse train $w[n]$, in the figure-b, by its (DFS) discrete-Fourier series sum:
$$ w[n] = \frac{1}{M} \sum_{k=0}^{M-1} e^{j \frac{ 2\pi k}{M}n }$$
If you cannot see why this is so, consult a DSP textbook for DFS.
In the third step observe that you can treat $v[n]$ as an expanded version of $y[n]$. Based on this, express Z-transforms of the expanded signal $v[n]$ in terms of Z-transform of $y[n]$ which is a lot easier to derive than the inverse operation of downsampling;  more specifically:
$$V(z) = Y(z^M)  \longrightarrow  Y(z) = V(z^{\frac{1}{M}})$$
In the fourth step, we shall express Z-transform $V(z)$ of $v[n]$ in terms of Z-transform $X(z)$ of $x[n]$, based on the following two properties of Z-transforms:
$$ \sum_k x_k[n] \longrightarrow \sum_k X_k(z)  ~~~ , ~~~ \text{linearity}$$
$$  e^{j\omega_0 n} x[n] \longrightarrow X(e^{-j \omega_0} z)  ~~~ , ~~~ \text{complex modulation}$$
Then express $v[n] = x[n] w[n]$ and apply the above two properties:
$$ v[n] = x[n]w[n] = x[n] \frac{1}{M} \sum_{k=0}^{M-1} e^{j \frac{ 2\pi k}{M}n } = \frac{1}{M} \sum_{k=0}^{M-1} e^{j \frac{ 2\pi k}{M}n } x[n] $$
$$ V(z) = \frac{1}{M} \sum_{k=0}^{M-1} X( e^{-j \frac{ 2\pi k}{M}} z) $$
In the final step, merge the above result with the result of step-3 by replacing $z$ with $e^{j\omega}$ to convert the result from Z-transform to discrete-time Fourier transform:
$$ Y(z) = V(z^{1/M})$$
$$ V(z) = \frac{1}{M} \sum_{k=0}^{M-1} X( e^{-j \frac{ 2\pi k}{M}} z) $$
$$ Y(z) = \frac{1}{M} \sum_{k=0}^{M-1} X( e^{-j \frac{ 2\pi k}{M}} z^{1/M}) $$
$$ Y(z) = \frac{1}{M} \sum_{k=0}^{M-1} X( e^{-j \frac{ 2\pi k}{M}} e^{j \omega/M}) $$
which is the relation between the discrete time Fourier transforms of the input $x[n]$ of a compressor and its output $y[n] = x[Mn]$ as:
$$ Y(e^{j\omega}) = \frac{1}{M} \sum_{k=0}^{M-1} X( e^{j \frac{w- 2\pi k}{M}}) $$
Interpretation of the final result is as follows: The discrete time Fourier transform $Y(e^{j\omega})$ of the the downsampled signal $y[n] = x[Mn]$ is given by the shifted and scaled sum of $M$ copies of DTFT $X(e^{j\omega})$ of $x[n]$ scaled in magnitude by $1/M$ and expanded in frequency by $M$ and each copy shifted by $2\pi$ to the right.
This expansion in frequency $\omega$ by M of the DTFT $X(e^{j\omega})$ can be pictured as follows: Consider the original input frequency spectrum $X(e^{j\omega})$, which is bandlimited to $\omega = \pi/M$ for preventing aliasing due to subsequent downsampling. Being a DTFT, by default, It's periodic by $2\pi$ having its centers at $\omega = ... -2\pi, 0 , 2\pi, 4\pi, 6\pi ...$. Now expanding the frequency axis by M yields a function $X(e^{j\omega/M})$ such that it's periodic by $2 M \pi$ and it has its centers at $\omega = ..., -2M\pi, 0 , 2M\pi, 4M\pi, 6M\pi ...$. This spectrum has complete zeros at those frquencies $\omega = 2\pi, 4\pi, (M-1)2\pi ...$, due to frequency axis scaling, and those M copies of these expanded spectrums are added to each other after being shifted by $2\pi k$ to the right, for $k=0,1,2,M-1$, and therefore filling the gaps in between $\omega = 0$ and $\omega = M2\pi$ of the spectrum of $X(e^{j\omega/M})$, hence having a period of $2\pi$, $Y(e^{j\omega})$ is now a proper DTFT.


Deriving vector paths from pen strokes I want to calculate vector paths that approximate pen strokes from handwritten documents. My tool needs to run in batch mode, unlike hand tracing or Illustrator's "Image Trace" feature.

The potrace vectorizer will only take black and white input -- but applying a color threshold to the scans destroys the useful detail in their smooth, blended edges.

I had some success with the SUSAN edge detector. It can take greyscale input and its output can be reduced to black and white. Here's a liberal set of edges found around the perimeter of the word "intention":

Now, how can I extrapolate the pen's possible path across the semi-enclosed area? I am aware that there will be ambiguities; potrace has several strategies for resolving them.

More broadly, are there better approaches to this problem?

...how can I extrapolate the pen's possible path across the semi-enclosed area?

There are infinitely many ways that this path could have been generated. Therefore, you will either have to make assumptions or use extraneous information for something that resembles "...the pen's possible path..." (which would not be too far away from the way potrace does it).

In any case, instead of applying edge detection, I would suggest that you apply skeletonisation first to get rid of the stroke's width and then try to work out the "evolution" of a much simpler path. The stroke's width generates what looks like a closed path which is incredibly hard to work out because of ambiguities (do I trace this clockwise or counter clockwise?)

Even when the skeletonised output is given, you still need extraneous information.

For example:

Where did the pen land on the paper for the first time?

You could select the "left-most" pixel for that but cursive writing can result in heavily slanted letters where this assumption is false. Checkout the "y" on the word "yesterday" in this image for example (image taken from the Wikipedia article linked earlier).

"Which process generated this curve?"

This is basically what you are looking for and possibly cursive writing could be expressed as some form of minimisation problem. That is: Find the shortest length curve that links as many of the given points as it is possible. Notice here, "curve", singular, one curve.
An algorithm that does this is Alpha Star. Alpha star operates over a graph data structure and finds a path that is optimal given a heuristic. Usually, this heuristic is an expression of "cost". Say for instance, the "cost" of establishing an edge between two points is proportional to their distance (definitely) but also how close that edge would bring us to the final destination. The algorithm basically calculates this for all transitions and then works out the path with the highest combined heuristics.
In your case, the graph data structure are all the points that result from skeletonisation connected at least with the very simple rule of shortest distance. So, you start at your "starting point", scan its neighbourhood for the nearest ink point, connect it, jump to the nearest and repeat. More than one pixels could be "nearest" in this case, which is fine at this stage, you connect them and then let A* tell you how to resolve a "fork".
In your case, the cost is obvious: The distance between pixels. But the heuristic is possibly where the game is won or lost. You can use the "how close it brings us to the end of the curve" but this means that you have an idea of where the curve ends.... Which brings us to even more assumptions:

"Where did the pen lift from the paper?"

You can possibly identify far enough sequential points and mark them as the "end points". This would give you rudimentary start-end points for tracing the curve but where does this leave the dot above the "i"? Was it added while writing the "i" or at the end of a given "curve"?

For these reasons, I do not think that you would be able to recover the pen's movement without any extraneous data.

It would help even if you could make an assumption on the pen nib that the writing was made with because if it was a flat nib and it was held at an angle, it can give you information about its "orientation" and it is safe to assume that the writer was not rotating the pen in their hands. This gives you an additional criterion to that of the nearest distance. Nearest distance and smooth angle transition.

A more accurate method would be to record cursive writing movements for each letter and then recognise each (visually) and "fit" the movements to the point data with "transition curves" between the letters. This would give you the possible way a given writer could have written a particular script. You can capture this data with lots of different devices (tablets, e-notebooks).

Hope this helps.


Z-domain transfer function to difference equation So I have a transfer function $ H(Z) = \frac{Y(z)}{X(z)} = \frac{1 + z^{-1}}{2(1-z^{-1})}$. I need to write the difference equation of this transfer function so I can implement the filter in terms of LSI components. I think this is an IIR filter hence why I am struggling because I usually only deal with FIR filters. I have tried to simplify the filter, and I get:

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

This gives me the gain (K = 0.5) and the poles as +1, and the zeroes as -1, hence filter is stable.

Can anyone help me? I usually divide through by the denominator, and hence get the difference equation, but I can't in this case. Then it's just a  case of looking at the difference equation and implementing the filter with delays, multiples, adders etc.

$$\begin{align*}\dfrac{Y(z)}{X(z)} &= \dfrac{1+z^{-1}}{2(1-z^{-1})}\\
\\
2(1-z^{-1})Y(z)&=(1+z^{-1})X(z)\\
\\
Y(z) -Y(z)z^{-1}&= \frac{1}{2}X(z) +\frac{1}{2}X(z)z^{-1}\\
\\
y[n]-y[n-1]&=\frac{1}{2}x[n] + \frac{1}{2}x[n-1]\\
\\
y[n]&=y[n-1]+\frac{1}{2}x[n] + \frac{1}{2}x[n-1]\\
\end{align*}$$

You'll need at least one register for the feedback of the output (the y[n-1] term).

You might as well call that register the accumulator, as this is essentially an integrator, integrating up the average of the two most recent input values.  Integrators are not stable, so you may want to clip the output at some point instead of letting it roll over.


Which sections to cover for my graduation thesis on blind source separation? I'm not sure if I'm on right place to ask a simple suggestion, so please be kind and don't downvote my question. I'm not asking someone to do it for me, I just need short guidance how should I proceed.

Currently, I'm writing a study of Blind source separation using eigenvalue decomposition. Since I'm not good at writing, I have done so far the needed programming part in Matlab. The study must be composed of ~50 pages theory about the problem. Theory should be in terms of explaining the BSS, BSS methods, ICA, and EVD which is the main topic. Can you suggest what to include in my writing since BSS has a lot of methods and it is difficult for me to decide and follow what to write and what not write. Can someone outline the order of things I should write and mention in my study?

I'm stuck at what should I write next, so far I have:

Introduction to BSS
Definition of Independence of signals
Independence and correlation
PCA and other prewhitening methods:
        1. centering
        2. whitening

I'm using books like: Independent Component Analysis - A Tutorial Introduction Handbook of Blind Source Separation Independent Component Analysis and Applications 2010 - Pierre Comon, Christian Jutten

Thank you,

I've taken this course.
I think using these two papers plus those books you mentioned would be so useful in writing a report, thesis , ... .

Blind signal separation: statistical principles

Independent Component Analysis (A Tutorial)

Also you could download my professor's thesis which is available online in this page


How do i select my frequency range to plot when fft is made on time domain I have a vibration signal that i need to convert from time domain to frequency domain using fft in python.

from scipy.fftpack import fft
yf =  fft(df["x"])
plt.plot(df["x"])

And i would like to plot it without DC value at 0Hz. So i neglected yf[0] and took N/2 frequencies to plot as per Nyquist theorem.

But I am having little trouble how to calculate x axis frequencies. how do i select my range of frequencies to plot? Also the amplitude is plotted from 0 to 25. Does this amplitude (Frequency domain) corresponds to sum of amplitude of a particular signal in time domain? How do i know which amplitude from time domain corresponds to amplitude in frequency domain?

The Discrete Fourier Transform takes a signal from the discrete time domain to a discrete frequency one. In the process, it is possible to evaluate the DFT in less (or more) "frequency bins".

In the time domain, the signal is indexed by a variable that denotes discrete samples. For example $x[n]$ where $n$ relates to physical time by $n \cdot T_s$ where $T_s$ is the sampling period.

In the frequency domain, the spectrum of the signal is indexed by a variable that denotes discrete frequency. For example $X[k]$ where $k$ relates to physical frequency by $k \cdot \frac{Fs}{N_F}$ where $N_F$ is the length of the DFT.

So, given a signal that was acquired at $Fs = 44100$ Hz in a buffer that was 1024 samples long (that is the $x[n]$) which was then transformed by a 512 point DFT (that is the $X[k]$), the correspondence between $k \in [0 .. 511]$ and physical frequency $f$ would be $f_k = k \cdot \frac{44100}{512}$ or a spacing of 86 Hz.

You can then use this correspondence to choose which part of the spectrum you want to plot.

As a side note, your accelerometer board is pausing data logging to write the acquired buffer to an SD card or send it over some interface elsewhere. This creates very short zero value gaps, which, overall, are going to distort your spectrum into something that looks "wavy" in a very distinct way. In fact, your spectrum is effectively the spectrum of your vibration data as modulated at 100% by a square wave at whatever frequency is implied by the length of the buffer and the speed of the SD Card / interface of your data logger. Depending on how short your recordings are going to be and how they are timed, you are going to have to do something about this, otherwise your spectrum will not be the "real" one.

Hope this helps.


Uneven background subtraction: Rolling ball vs Disk tophat A classical technique in still images (e.g. fluorescence microscopy images) to remove uneven illumination and isolate bright blobs is to use morphological operation such as the top-hat transform.

For instance the Rolling-ball algorithm[1] uses a ball as a structuring element and performs the top-hat transform[2].

The structuring element is a grayscale element that extends the disk structuring element in 3D. For instance in MATLAB:

se2 = strel('ball', 3, 3);
     0    0.4999    0.9997    1.4996    0.9997    0.4999         0
0.4999    0.9997    1.4996    1.9995    1.4996    0.9997    0.4999
0.9997    1.4996    1.9995    2.4994    1.9995    1.4996    0.9997
1.4996    1.9995    2.4994    2.9992    2.4994    1.9995    1.4996
0.9997    1.4996    1.9995    2.4994    1.9995    1.4996    0.9997
0.4999    0.9997    1.4996    1.9995    1.4996    0.9997    0.4999
     0    0.4999    0.9997    1.4996    0.9997    0.4999         0

My question is: Is it so much better than simply using a binary disk as a structuring element? What are the use cases where a rolling ball is way better, or even required?

My experiments with comparing the outputs of the two do not show a striking difference between the two, but the literature behind the Rolling ball is very abundant.

[1] Stanley Sternberg, "Biomedical Image Processing", IEEE Computer, January 1983

[2] Here is the source code for a Java implementation in ImageJ: http://rsbweb.nih.gov/ij/plugins/download/Rolling_Ball_Background.java

I guess the basic question here becomes - what difference does non-flat structuring elements make w.r.t flat structuring elements ?

From the definition of dilation one can see that the structuring element can define a support (like in flat structuring elements) or another function itself (gray scale structuring elements). Look at wiki which defines the dilation operations. The grayscale or numerical version performs a Minkowski Sum. The dilation operation here basically looks for the max on the sum of the two functions. With grayscale dimension one is basically looking at the topographical features of the image level. Its called the rolling ball algorithm since it shaves off non smooth parts of the image which are less convex that the sphere structuring element one defines. The top hat is non flat and basically removes features in each level set and has no effect on the topographical height directly. It depends on the image and problem how these 2 operations are used. One can also write structuring elements to look for connected maximas or minimas(Jump connection), smooth regions(Smooth Connection) refer.


Audio phase shift using Hartley Transform I want to shift the response phase of an audio signal (float array, length N, value range -1..1). I'm about using Fast Hartley Transform to transform it to the frequency-domain, process phase shift, then convert back to the time-domain and output the result.
I'm now can do the Hartley transform, and I also can calculate the phase of kth bin:

real(k) = (x(k) + x(-k))/2
imag(k) = (x(k) - x(-k))/2
phase(k) = atan2(imag(k), real(k))

I want to shift the phase(k) by an α degree. What should I do?

I think your notation is a bit fishy but I'll use it. Assuming that your $\text{phase}(k)$ calculation is correct, do a complex multiplication of $\text{real}(k) + i\ \text{imag}(k)$ by $e^{i\alpha} = \cos(alpha) + i\sin(alpha)$. That can be done using four real multiplications and a few add/sub. Then you need to reverse the process of how you obtained $\text{real}(k)$ and $\text{imag}(k)$:

$$x(k) = \text{real}(k) + \text{imag}(k)$$
$$x(-k) = \text{real}(k) - \text{imag}(k)$$

Finally retransform.


How does Signal Analyzer measure EVM I understand that EVM is a root mean square of error vectors. If I have transmitted symbols and received symbols (once equalized), taking a root mean square would give me EVM. I believe this is actual EVM (we are using actual tx symbols and mapping them to equalized received symbols and then computing EVM).
How does a signal analyzer (which doesn't know tx symbols) computes EVM?
Does it calculate the EVM  of each equalized received symbol with a closest ideal symbol on the constellation and then average it out to give final EVM value?

Does it calculate the EVM of each equalized received symbol with a closest ideal symbol on the constellation and then average it out to give final EVM value?

You very rarely test EVM on a low-SNR signal that's not according to any standard, especially because if that signal doesn't adhere to any standard, there would be very little for your signal analyzer to do equalization with.
If you buy a signal analyzer, you can often buy extensions that allow you to estimate properties of some specific communication / broadcasting standards (LTE, DVB-T, VDSL,…). Such standards have channel coding, with which you can (most of the time) recover incorrectly received symbols, and then actually calculate the error between reception and actually received symbol.
In other words, to know what your signal analyzer does in your current situation, you'll have to read its documentation.


Can I control per-pixel contrast in an image? Is there any way to control per-pixel contrast? E.g., a way to set an upper or lower limit on the contrast between pixels?

The closest thing I've found is a high-pass filter of an n-pixel radius, inverted and applied as an overlay layer in Photoshop. Other kinds of filters, such as median, blurs, and traditional contrast adjustment, produce a result which is similar to my goal in some circumstances, but for the whole image, rather than only pixel pairs over a given contrast threshold, or proportionally based on contrast.

Does this concept have a name? Is there an easier way to do what I'm attempting?

[Edit: And is there any free or affordable software which can be coerced to do this?]

By per-pixel I assume that you mean only the surrounding pixels should influence the pixel value. Unfortunately changing the pixel value in response to one adjacent pixel means the contrast with the others could also change. Kind of like a 'chain reaction' of changes.

For example, saw we had $I = [0,7,11]$ and we wanted a maximum contrast of 5. Going from the left we change 7 to 5, which then requires 11 to be changed to 10.

Maximum Contrast

One way per pixel maximum contrast can be done is by using morphological operators. Specifically, performing a dilation (or erosion) with a pyramid structuring element based on the maximum contrast value, $c$. For dilation, the structuring element is

... ... ... ... ... ... ... ... ...
... -3c -3c -3c -3c -3c -3c -3c ...
... -3c -2c -2c -2c -2c -2c -3c ...
... -3c -2c -1c -1c -1c -2c -3c ...
... -3c -2c -1c   0 -1c -2c -3c ...
... -3c -2c -1c -1c -1c -2c -3c ...
... -3c -2c -2c -2c -2c -2c -3c ...
... -3c -3c -3c -3c -3c -3c -3c ...
... ... ... ... ... ... ... ... ...

The size of the structuring element will depend on $c$ and the range of pixel values. Assuming pixel values are in [0,255], size will be N x N pixels where N = 2*floor(255/c) + 1, e.g.

c = 100

-200 -200 -200 -200 -200
-200 -100 -100 -100 -200
-200 -100    0 -100 -200
-200 -100 -100 -100 -200
-200 -200 -200 -200 -200

Dilation increases pixel values to enforce the maximum contrast, erosion decreases pixel values. In out example above, dilation would give $I = [2,7,11]$, erosion would give $I = [0,5,10]$.

Minimum Contrast

Enforcing a minimum contrast is not going to give you good results. Imagine having a minimum contrast in a flat area, it will turn it into a ramp. For this operation it is better to look at filters that enhance edges etc.

MATLAB code

MATLAB code to do this:

clear all;
    close all;

% Test image
I = randi(256,10,10)-1;
imagesc(I); colorbar; title('test image'); pause;

% Constrast settings
c = 10;
maxPixelValue = 255;

% Structuring element
halfN = floor(maxPixelValue/c);
N = 2*halfN + 1;
[x,y] = ndgrid(-halfN:halfN);
stDH(:,:,1) = abs(x) * -c;
stDH(:,:,2) = abs(y) * -c;
stDH = min(stDH,[],3);
stD = strel('arbitrary',ones(N),stDH);
stE = strel('arbitrary',ones(N),stDH);
imagesc(stDH); colorbar; title('structuring element'); pause;

% Dilation result
Id_max = imdilate(I,stD);
imagesc(Id_max); colorbar; title('max contrast (dilation)'); pause;
imagesc(I-Id_max); colorbar; title('diff with original'); pause;

% Erosion result
Ie_max = imerode(I,stE);
imagesc(Ie_max); colorbar; title('max contrast (erosion)'); pause;
imagesc(I-Ie_max); colorbar; title('diff with original'); pause;

%% Min contrast
stD = strel('arbitrary',ones(N),-stDH);
stE = strel('arbitrary',ones(N),-stDH);

% Dilation result
Id_max = imdilate(I,stD);
imagesc(Id_max); colorbar; title('min contrast (dilation)'); pause;
imagesc(I-Id_max); colorbar; title('diff with original'); pause;

% Erosion result
Ie_max = imerode(I,stE);
imagesc(Ie_max); colorbar; title('min contrast (erosion)'); pause;
imagesc(I-Ie_max); colorbar; title('diff with original'); pause;


Need critical help: How to detect and distinguish two very similar looking signals? Hi guys I have a really tough signal processing question here.

How do you detect and distinguish two very very similar looking waves?

I need to distinguish between these two signals for an online application and the techniques so far have all failed. These two signals have very similar and often unpredictable amplitudes. Almost identical frequency. So the only way I can attempt to distinguish between the two is using time domain features.

Signal 1:

Signal 2:

The data is streamed live in every single second (at 128 samples per second). My approach is to:

use a window of size 128 samples to obtain the data
Calculate the time domain features of the data
features such as:

derivative between the mean of the current 128 samples and previous
128 samples,
derivative between the max point of the current 128 samples and
previous 128 samples
skewness of the data
kurtosis of the data
deviation from the mean

If derivative goes down, then we say it is signal 1 and if derivative does down we say it is signal 2

But the problem is that the two signals look very very similar. So if the window is incorrectly placed (i.e. on the rising or falling phase after the initial signal), then it can give me the complete opposite result!

Can someone give me some suggestions as to what features I should use or what technique I should use to distinguish between this two very similar looking signals!

As Conrad pointed out, a correlator is probably your best bet.

The correlation of a signal with itself (also known as its self-similarity) is larger than its correlation with any other signal (except for a constant factor related to the signals' energy).

In your case, you would implement two correlators, one for Signal 1 and one for Signal 2. Then, you'd look for one of the correlators' output to be larger than a certain threshold. You may have to run a few tests to determine the threshold.

One nice benefit of using correlation is that you don't need to synchronize your time window to the signals. Just be sure to accumulate the correlation calculations from one time window to the next.


Adaptive filter to scale and phase shift two sensors output One way of separating downgoing and upgoing wavefields in offshore seimic processing is to add signals from hydrophone and vertical component of the geophone (they are co-located). Hydrophone only registers a change in the pressure whereas geohphone as well as registering a change in seismic field also reacts to the direction of arriving wave, ie. if it is coming from above or below (sensors are deployed on sea floor). But the scale of the two are different as they record two different physical quantities. In addition, there is a phase difference between the two. To do the summation, I first need to re-scale one to the other and phase shift. The phase shift is slightly frequency dependent. I would like to make a filter that does the above on a period of data and then use the filter on all incoming data. The figure below shows an example of an upcoming signal (hydrophone signal has been downscaled for the exmple. In reality it is of much higher amplitude). If the two were of same scale and in phase then summing them will double the signal's amplitude while reducing the noise coming from above. I have some solution in time domain (assuming constant phase shift and scaling factor) and frequency domain by multiplying one's spectra with the transfer function between the two (found from spectral division). I have also used 'nlms' adaptive filtering. But some noise still remains, enough to make trouble later on.
Wonder if anyone has a suggestion on how to make the filter or some Matlab code.
Thanks

It seems that you want to equalize one time series to the other or vice versa or some combination. The general problem is then (in shorthand)
$$
h_1{\Large *}y_1 \approx h_2{\Large *}y_2
$$
where $\Large *$ denotes convolution, $h_i$ are filters, and $y_i$ are your time series. If you simply equalize one time series to the other, this is
$$
y_1 \approx h_2{\Large *} y_2
$$

In matrix form, this becomes
$$
\left[\begin{array}{c} y_1[p]\\y_1[p+1]\\\vdots\\y_1[p+N-1]\end{array}\right] \approx
\left[\begin{array}{cccc}y_2[2p] & y_2[2p-1] & \cdots & y_2[0]\\
y_2[2p+1] & y_2[2p] & \cdots & y_2[1]\\
\vdots & \ddots & \ddots & \vdots \\
y_2[2p+N-1] & y_2[2p+N-2] & \cdots & y_2[N-1]
\end{array}\right]
\left[\begin{array}{c}h_2[0]\\h_2[1]\\\vdots\\h_2[2p] \end{array}\right]
$$
or more concisely
$$
\underline{y}_1 \approx {\bf Y_2}\underline{h}_2
$$
The least squares solution for the equalizing filter coefficient vector is
$$
\underline{h}_{2,LS} = \left({\bf Y_2^HY_2}\right)^{-1}{\bf Y_2^H}\underline{y}_1
$$

You can alternatively use LMS or RLS to do recursive updates to find/track $\underline{h}_2$.

One additional comment on equalizing one time series to the other: by choosing to equalize $\underline{y}_2$ to $\underline{y}_1$, we choose one spectral weighting - implicitly. If we instead equalized $\underline{y}_1$ to $\underline{y}_2$, we would have a slightly different spectral weight for the overall least squares fit. If the filters are relatively minor in magnitude variation over the signal band, the spectral weight difference is similarly minor, but if the equalizer filter response is quite different for the two cases (can just as well solve and analyze both), you will need to assess which is preferred/better and why.


Simple analog anti aliasing filter good enough if signal frequency is far from sampling frequency ? What about equivalent time sampling? I have one and only one frequency signal, a sine signal at f = 10 kHz that is sampled at Fs = 320 kHz.

I have read about applying a simple analog anti aliasing (AAA) filter then oversample, do some digital anti aliasing (DAA) and decimate to get the correct signal frequency.

I understand how this is useful if signals of interest are close to Fs/2 as it might be hard to find an AAA filter that has sharp enough roll-off, low cost, simplicity etc. And therefore DAA filter in combination with a simple AAA filter is often preferred.

But in my case I don't have a signal close to Fs/2 at all. So I was thinking maybe I can get a slower/cheaper ADC and do no oversampling, since I might be able to properly attenuate frequencies that might alias even using a simple AAA? Is what I am thinking correct?

Additional related question:
How would this work if I was to apply equivalent time sampling (ETS) technique? This means that the signal is 10 kHz and the sampling frequency appears to be 320 kHz as well, but 32 samples are taken over 32 periods of the 10 kHz signal to produce 1 period. So the sampling rate is actually 10 kHz + (10 kHz / 32) = 10.3125 kHz. What I am not able to wrap my head around: will Fs = 10.3125 kHz or 320 kHz? Will the same filter work for both cases? I assume Fs = 10.3125 kHz thus giving a problem as the transition area will be relatively very small: Fs - f = 3125 Hz.

Other thoughts:
I am using a DDS to generate the sine signal, it should be quite pure and maybe AAA is not "needed" / very important?

The sampling theorem states that you need at least "two samples per one Hz of bandwidth". Anti-aliasing filters are required when the bandwidth of the original analog signal is not well defined and a finite bandwidth needs to be enforced.

In your application that doesn't seem to be the case. If all you want to sample is really just an unmodulated steady-state sine wave, the bandwidth or your original signal is already limited. No anti aliasing filter is needed and you can simply sample at any rate slightly larger than twice the sine wave frequency. 22.05 kHz or 24 kHz in your case would work fine.

ETS is a bit tricky since you need some sort of a trigger mechanism and then a variable delay off that trigger. The equivalent sample rate depends on your delay steps and the whole process depends on the signal being periodic with high precision.

You actually CAN sample at a much lower rate. The theorem requires 2 samples per Hz of bandwidth. For a pures sine wave the actual bandwidth is almost zero so you can sample with a very low sample rate and still don't have any loss of information. You just need to manage the mirror spectra properly, i.e. your application needs to be aware of a frequency offset.


Why do we choose cutoff frequency to be taken for "1/sqrt(2) magnitude reduce"? Why don’t we instead choose (calculate) our bandwidth to be between two purple “dot’s”? Why don’t those purple dot’s represent points for “cutoff frequency”? Looking at this picture “frequency components” after the right purple dot (or before the left purple dot) would be attenuated (yellow lines + after cutoff frequency). I presume that red line I drew represents that components of a signal would neither increase or decrease in magnitude. Following that logic, if I want to send signal and (theoretically) receive same signal through this channel (having this response) than it’s better for signal to have only components between those two purple dots. Any components of signal after right purple dot or left purple dot would be attenuated and hence we won’t receive same signal (signal with same shape)?

I presume that red line I drew represents that components of a signal
  would neither increase or decrease in magnitude.

That would be a wrong presumption. The magnitude of the frequency response is (almost) never truly flat, and you would see this the more you zoom in. That's why you need to pick a convention of what exactly the "band edge" is.

As with most conventions, it's somewhat arbitrary. Half power has the advantage that if adjacent filters have the same cutoff frequency the total power will be flat. So if you add the output of an 1 kHz high pass and a 1 kHz low pass, the energy sum of the outputs will be flat with no dips or peaks in the cross over region.

-3dB is also not the only convention either: a Linkwitz Riley filter typically uses -6dB at the nominal cutoff. -3dB also doesn't make a lot of sense for parametric EQ filters and for high and low shelve filter.

So it all depends on context and one should never automatically assume -3dB.


I think I might have finally solved it myself.

Afaik, while it is true that the actual allpass filters are 1st order, the trick is to see that the original formula actually just says that the two chains of allpasses run interleaved. That is the reason for the z^-2 and the reason for the extra multiply by z^-1.
Hence it is perfectly reasonable to simply use z^-1 filters in the implementation, and feed the allpasses the data interleaved.

However, for plotting the response, I had to first compute the cascaded filter coefficients, and then finally apply the z^-1 delay to the A1 chain.

This is the result:

Which is exactly what I wanted.
This is the final code:

c = [0.036681502163648017, 0.13654762463195794, 0.27463175937945444, 0.42313861743656711, 0.56109869787919531, 0.67754004997416184, 0.76974183386322703, 0.83988962484963892, 0.89226081800387902, 0.9315419599631839, 0.96209454837808417, 0.98781637073289585 ];
K = length(c);
i=0;

b0=[1];
a0=[1];
b1=[1];
a1=[1];

while (i<K)

  beta0 = c(i+1);
  b0=conv(b0,[beta0,0,1]);
  a0=conv(a0,[1,0,beta0]);

  beta1= c(i+2);
  b1=conv(b1,[beta1,0,1]);
  a1=conv(a1,[1,0,beta1]);

  b1z1 = [0,b1];

  [A0,w] = freqz(b0,a0,512,2);
  [A1,w] = freqz(b1z1,a1,512,2);

  A=(A0+A1)/2;

  plot(w, 20*log10(abs(A))); % amplitude plot in decibel
  hold on;
  ylabel("Gain (dB)"), axis([0 1 -120 20])
  xlabel("Normalized frequency");
  grid on

  i+=2;
end

I think this should be the accepted answer, but if anyone has anything to add, please do.


Digital Downconversion of Signal DSP-Newbie here... Why are we using the ratio of the center frequency of the signal to the sampling frequency of the signal spectrum when we downconvert the signal?

I have an Spectrum with Carriers which reaches from $-f_{s}/2 ... f_{s}/2$ ($f_{s}$ is the Sampling-frequency) (see picture)

The downconversion/shifting of the spectrum is

$$e^{-j2\pi(LO/f_s)}$$

For a centre frequency of $F_C$ for a modulated signal $v(t)$.
$$
r(t) = v_I(t) \sqrt{2} \cos 2\pi F_C t - v_Q(t) \sqrt{2}\sin 2\pi F_C t
$$

Assuming a direct RF sampling at a rate of $F_S=1/T_S$ before downconversion, we plug $t=nT_S$ above.
$$
r(nT_S) = v_I(nT_S) \sqrt{2} \cos 2\pi F_C nT_S - v_Q(nT_S) \sqrt{2}\sin 2\pi F_C nT_S
$$

The expression $F_C nT_S$ above becomes $\frac{F_C}{F_S}n$ to yield
$$
r(nT_S) = v_I(nT_S) \sqrt{2} \cos 2\pi \frac{F_C
}{F_S} n - v_Q(nT_S) \sqrt{2}\sin 2\pi \frac{F_C
}{F_S} n
$$

Perhaps this is the ratio of the centre frequency to the sample rate you are asking about.

EDIT: Downconversion can be embedded in a downsampling operation as follows. In a conventional operation, we will have an LO downconverting the signal to baseband, filter it and then downsample it according to the symbol rate, as shown below.

Here, the filter is operating at a higher rate which is unnecessary when we have to throw $M-1$ out of every $M$ samples anyway. This is a hint for what we need to do.

First, we interchange the operations of downconvesion and filtering so our filter needs to be implemented in passband now. If the lowpass filter impulse response is $h[n]$, the passband filter response is
$$ g[n] = h[n] e^{j2\pi \frac{F_C}{F_S} n}$$

The interchange of downconversion and filtering is shown below.

Next, we slide the downsampling operation past the LO whose frequency gets multiplied by $M$, i.e., $F_C/F_S \cdot M$. We assume that
$$F_C = \frac{F_S}{M}$$
then our LO frequency after sliding it past the downsampler becomes
$$\frac{F_C}{F_S}M = \frac{F_S}{M}\frac{M}{F_S}=1$$
which implies that no downconversion is actually required!

However, remember that filter is still operating it at a higher rate. So we slide the downsampler past the filter now and implement it as a polyphase filter. The original spectrum is shown below.

This operation changes the spectrum such that all spectral replicas from $-F_S/2$ to $F_S/2$ end up at baseband overlapped on each other. The interesting part is that at the input of the polyphase arms, they possess a unique phase profile across the whole $360^{\circ}$ around the frequency IQ spectrum. A rough sketch is below.

This is only true in an exact manner for the sinusoidal signals, so the rest of the task for a signal with bandwidth is performed by the filter. When implemented as a polyphase filter, this filter operates with the delays matching those of the spectral replicas shown above. A final summation at the end cancels everything out leaving the desired spectrum.

This also implies that if the desired spectrum is at $F_S/M$ and we want to employ the initial lowpass filter $h[n]$ instead of the bandpass filter $g[n]$, a corresponding rotation is required after each polyphase arm to line everything up before the final cancelation.


How to Extract a Period of a Periodic Pulsed Signal? I rercorded a sequence of zeros and ones which is repeated many times. However, the first and the last repetitions may not be complete because I have to start recording at a random time, so the sequence can be started when I hit record and the sequence may have not finished when I stop recording.

I'd like to know how to extract, for instance, the 3rd period of the complete recording, because this is something I know it would have the complete sequence.

here is a example of the signal I'm working with. I computed the FFT and I think the frequency is 96000Hz (half of the sampling rate), but still I'm not sure how to extract a complete period from it.

How could I do this?

Natalia Molinero Mingorance, Welcome to the DSP community.

What you have is basically a shifted periodic signal.
Why? Because what you have is equivalent (Given many samples) of having a periodic signal which you started sampling at $ 0 $ and then we shit it to the real place you started sampling it.

Yea, there is the edge case, since it might not have  its end and beginning to make a complete cycle.

But let's assume our model is perfect, what would you see in Frequency Domain (DFT)?

Well, since shifting in time only changes the phase it means you won't see any difference.

So how would you extract the period of a signal you started sampling at the beginning of a period in the Frequency Domain?
Yep, just looking at the Frequency in the Frequency Domain.

If you post the files of the samples we'll be able to to see this together.


Envelope following using FFT bin power -> Low pass effect on the modulation transfer function I am deriving the temporal envelope of a signal (in a particular frequency band) by summing some FFT bin powers and taking the average, I do the FFT frame by frame, windowing and some overlap.

Please help me understand the following: when the envelope of the signal changes slowly (say 5 Hz), I can track the envelope fine. However when the envelope is modulated at a higher frequency, say at 100 Hz, the dynamic range and amplitude of the detected envelope are very reduced.

So the Modulation Transfer Function (MTF) of this method has a low-pass look.

I am trying to understand the reasons for this low-pass roll-off effect and the maths behind it. I see what's happening but I can't clearly conceptualise the underlying cause . It should also be possible to come up with a number regarding the slope of the MTF.

Thanks in advance!

GP,

Edit 1 : Thanks for your answer. I can't post images cause I'm a new user apparently.. so here is a link: Link . you see on the left the input signal is 4 kHz modulated by 5 Hz and I take some fft bin powers to track the 5 Hz envelope and it works fine. But at 75 Hz modulation (right), the amplitude and dynamic range of the envelope are reduced, and it gets worse as you increase the modulation freq (low pass effect). This is what I am trying to get my head around: where is this effect coming from ? Is it because I destroy phase information when I take the power of the bins ? I apply a hann window on the time domain input frames, before computing the FFT.

Any set of FFT framings of the data is a form of averaging filter or FIR filter, and thus will have a low-pass filtering effect over some envelope frequency range.  You can understand this effect by treating the FFT window as a FIR filter kernel and calculating its frequency response.

If the width of your FFT window is somewhere near the period of 5 Hz, then you will be in an area of filter roll-off.  You could try a shorter window, compensate for the window roll-off, perhaps using a "flat-top" window if your frequency of interest is well within the window's passband, or use another form of enveloper follower which does not have an intrinsic roll-off in the area of interest.


Determine Sound Intensity in SPL I am recording ultrasonic sound using "Earthworks M50 50kHz Measurement Mic" and the sound card is "Rolando Quad-Capture USB 2.0 Audio Interface". The recording was done using Matlab into a wav file. Can anybody explain how can we determine the intensity of sound. What does the aptitude of a wav file represent? Is there any formula that we can use to find out SPL using the above information?

Since you have a measurement microphone, it should come with a manual that contains all necessary information to characterise it (if not with all this information in digital format too).

One crucial piece of information is its Sensitivity which is expressed in V/Pa (Volts per Pascal), at a specific frequency. This figure allows you to convert VOLTAGE (which is what you are measuring at the output of the mic) to PRESSURE (which is what the transducer responds to).

The sensitivity of the microphone, in combination with its frequency response will tell you how much VOLTAGE does the microphone produce for a reference PRESSURE level at a given frequency.

So, this mic states 36mV/Pa and is relatively flat from 10Hz to 10kHz. This means that a sound wave of Sound Pressure Level (SPL) of 1 Pascal falling on the mic's transducer will cause a 36mV signal to be generated at its output at any frequency between 10Hz and 10kHz.

Outside of these ranges, the curve is slightly lowered which means that the sound wave's SPL will still be 1Pa, but the output of the mic will not be 36mV but 36mV times the factor that is obtained from the frequency response.

So, now we have a way of converting between SOUND PRESSURE LEVEL (in Pascal) and VOLTAGE. What about what is hapenning in the sound card?

The sound card will have a pre-amplifier / amplifier to amplify the signal from the microphone and this amplifier will have a frequency response on its own which will modulate the frequency response of the mic. In other words, because of its construction, the amplifier will amplify different frequencies by slightly different factors, just as the microphone did.

To cut a long story short: A sound card's LINE level is 1.736V. The amplifier of the sound card is set to 1 and you are recording with 16bit resolution. This means that the positive peak of 1.736V is converted to the digital value of $2^{15}=32768$.

What value does 36mV produces? $\frac{36 \cdot 10^-3 V}{1.736 V} \cdot 32768 = 680$

Right, so when the WAV file reads a sample of 680 at an amplification level of 1, then we know that the microphone was producing 36mV and through its sensitivity we can infer that 1Pa of sound pressure level was falling on it.

So, to convert between your sample value and the sound pressure level that caused this sample value you need:

The microphone's sensitivity
The microphone's frequency response
The sound card's LINE IN Spec
The sound card's amplification level
The wordlength that was used during sampling.

Hope this helps.


Why is anisotropic diffusion useful in image processing? I am working on anisotropic diffusion and the two coefficients proposed by Perona & Malik.

I wanted to know what is the use of diffusion in image processing? Why is anisotropic diffusion important and in which fields is it commonly used?

The anisotropic diffusion algorithm by Perona and Malik is the pioneering work in partial derivatives equations (PDE)-based denoising.

It applies the law of diffusion on pixel intensities to smooth textures in an image. A threshold function is used to prevent diffusion to happen across edges, and therefore it preserves edges in the image. (Unlike for instance gaussian blur filter.) This makes it very interesting if you want to remove noise, but do not want to smooth out the edges of your image, for instance if you want to use these edges to segment the image, without being perturbed by the noise.

Numerous efforts have been made that build on it, improving or extending it.

Now for where it is used, I only have a limited culture. I can cite two

Image analysis in the Life-Science field (where I work): the images that you can get out of a microscope are extremely noisy, and most of the time, it is by construction. Automated image analysis of these data often involve segmentation, for which you sometimes PDE-based algos.
Video games! For instance, try playing Mass Effect (the first one at least).


In MP it is very important which dictionary is used. For example, if the atoms are orthogonal The problem will converge faster.

Another interesting example, where two entries of the dictionary are the same, you can wither find both corresponding values in the sparse signal as 1 and zeros or as 0.5 and 0.5. You can find any number of representation for such an example. Which one will the algorithm converge to?

In general, the picking of the dictionary $D$ is a very important point in the MP algorithm. I recommend you take a look in the following link as well as researching a bit more yourself. The Fourier matrix, for example, is a common choice for the dictionary.

To answer your question, while I could not find a problem with your implementation of the algorithm, the data you are using might be problematic.


If two signals have different distributions that you know and can model, can you use that to separate them? For example, if I have a mixed signal composed of signal A added to signal B, and I know that the histogram of signal A is non-Gaussian in the time domain, and that the histogram of signal B is Gaussian in the time domain, is there a way to isolate either A or B using this information?
I understand that I can use MLE or some similar method to fit a Gaussian to signal B, but I'm not sure where I'd go from there. My initial idea was to create a Gaussian filter that matches the distribution of signal B, and use that to separate B from A+B, but I'm not sure how I would implement that in practice, or if that even makes sense.
Any insights would be appreciated! Thanks

If two signals have different distributions that you know and can model, can you use that to separate them?

Statistically, yes. Does that mean you can be certain at any point in time? No. But it allows you to transport data.

For example, if I have a mixed signal composed of signal A added to signal B, and I know that the histogram of signal A is non-Gaussian in the time domain, and that the histogram of signal B is Gaussian in the time domain, is there a way to isolate either A or B using this information?

Not to separate them, but to use both to transport data independetly: For example, you can switch the Gaussian source on and off, and detect the presence of it.

I understand that I can use MLE or some similar method to fit a Gaussian to signal B

Uff, you don't "fit" a Gaussian to something that has a distribution like a Gaussian. I know this sounds like nitpicking, but in these contexts, it's very important to make a distinction between the shape of the signal and the signal's amplitude distribution.

My initial idea was to create a Gaussian filter that matches the distribution of signal B, and use that to separate B from A+B, but I'm not sure how I would implement that in practice, or if that even makes sense.

Doesn't make sense, far as I can tell.
I'm afraid you'll have to look into information theory a bit, and there into what mutual information is, how we model that, and noise.
You can then make general statements on the ability of A to transport data in the presence of B, considering the latter "noise". It also works the other way around!
This directly leads to the successive cancellation decoder, where you first (try to) decide the better receivable variable, then subtract the decision from the observation, then decode the other. For example:
Your A is discretely uniformly distributed on ${-2, -1.9, -1.8, \ldots , -0.1, 0, +0.1,\ldots +2\$, your B is Gaussian with zero mean and variance 0.01. You observe  Y = A+B = -1.52; you decide that by far the most likely value for A is -1.5, and with that knowledge you can infer that B must have been -0.02.
You needed distribution knowledge of A and B here to decide which of the two to decode first, and then to make a maximum-likelihood decision on A.
This successive method gives you a lower bound for your ability to get data from A+B given their distributions. However, that's not the same as a joint decoding. General, for a maximum likelihood estimator, you'd need to know the joint distribution of A and B under A+B; then you can, for any A+B, tell the most likely A and B. There's to my knowledge no general formula for this decoding – it really depends on the distributions. Also note that if A and B are statistically independent, this joint density will fall apart into the product of the individual densities. Note, however, also that usually, both A and B will employ channel coding, so that these likelihoods are both time-variant and depend on other decisions – welcome to the wonderful world of coding theory.


Filling hatched areas in image I'm trying to extract / fill the walls in a floor plan, the walls are filled with a (cross-)hatch pattern at a more or less fixed angle.

I've managed to remove most of the "noise" by masking the DFT so I have both the +45 deg lines and -45deg lines in separate images.

The boundaries of the walls are quite obvious to the human observer but somehow I'm stuck at segementing those lines into groups, i.e connected components.

What I could think of is to apply a run-length encoding and then looking for subvectors of length 4 which fullfill a certain ratio of black and white. Or a cross-correlation of a pre-defined hatching pattern. The issue I have with these methods is that I'd have to introduce a parameter which has to be tuned for every floor plan, this is what I want to avoid.

So I'm asking you guys whether you could give me some input how to approach this problem.

As mentioned in another answer, image erosion is quite efficient to fill the gaps between the wall hatched lines. The key point is to choose a structure element that is big enough to allow connection to be done.

You should use a structure element with both horizontal and vertical dimensions of at least half the distance between hatched lines (try different values!). You may end with a final erosion with a small structure element to remove remaining noise (the thin and lonely lines in your case).

Then it is your responsibility to detect and refine the borders ;)

Here is attached the result I obtained with a square structure element of 19x19 pixels on a sample of the image you linked (logical OR one) assuming that the distance between lines is 19 pixels:

Hint: The output is smoother with a larger structure element but artefacts appear at the corners.


Choosing sampling frequency and base interval to eliminate leakage I'm studying for an exam and have some trouble with an exercise.
So I have a continuous-time signal that is

$$s(t)=\sin(2\pi f_{s1} t) + \sin(2\pi f_{s2}t)$$
with $f_{s1} = 1.4\text{ kHz}$ and $f_{s2} = 1.45\text{ kHz}$.

After sampling the time continuous signal I'll have to get the DFT of the discrete signal. The task is to chose the sampling frequency and the base interval of the discrete signal so that there will be no leakage-effect after calculating the DFT.

I also know that $\frac{N}{\text{period}}$ has to be an integer.

How would you approach something like this? I don't expect to be spoonfed the answer, hints would be nice though.

The key attribute you need to consider is for a DFT not to have leakage a sinusoidal signal must have a whole number of cycles in the sample frame.  Looking at your frequencies, clearly a one second interval satisfies this criteria, but it is not the shortest interval to do so.  The base interval is the shortest one that you are looking for that has a whole number of cycles for each signal.

Notice, this does not determine the sampling frequency.  Your $0th$ sample and your $Nth$ sample have to occur at the same relative point in your repeat pattern.  To get this, once you have found your base interval in continuous time, divide it into $N$ pieces to get your sampling frequency.  Any choice of $N>1$ will work to not have leakage, but if you want the bin indexes to represent the actual frequencies and not aliases you have to select $N$ large enough so that the higher frequency signal is at or below the Nyquist frequency.  Another way to say this is that you have to have at least two samples per cycle.

Followup:

This is the spoonfeeder version so the OP should stop reading here until they've tried to solve it using my clues themselves.

Honestly, I think Fat32's solution is a little overcomplicated, but impressive.

Let's call the two component sinusoidals Tone1 and Tone2.  Tone1 has 1400 cycles in a second and Tone2 has 1450 cycles in a second.  Observe:

$$ \frac{1400}{1450} = \frac{50 \cdot 28}{50 \cdot 29} =  \frac{28}{29} $$

From this we can deduce that it takes Tone1 28 cycles and Tone2 29 cycles to span the base interval of $1/50$ of a second.  Thus any DFT on this interval, no matter the sampling rate, will have non-zero bin values at bin 28 and bin 29 (or their aliased equivalent) and at the corresponding mirror image bins in the negative direction.

If $N \ge 29$ then it will be bins 28 and 29, otherwise you have to take modulo N to figure out the aliased bin index.  In order to be below the Nyquist limit $N \ge 58$.  Note that if either tone falls on the Nyquist bin (or the DC bin), the results may be zero or all the way up to the actual amplitude depending on where the sampling occurs within the sine wave.  For the OP's equation, the result will be zero since the sampling will occur at the zero crossings.  There still won't be any leakage.

For Fat32's first example, a samping rate of 4000Hz was used.  For an interval of $1/50$ of a second it means $ N = 80 $.  As per the first diagram, you can see that bins 28 and 29 are non-zero, and the mirror image bins 52(=80-28) and 51(=80-29) are also non-zero.

Fat32's second example is actually over an interval of 1/25(=15/375) of a second, so it is twice as long as the base interval.  Thus Tone1 has 56 cycles in the frame and Tone2 has 58.  Since N=15 in the example, aliasing does occur.  The bin indexes are then 11(=56 mod 15) and 13(=58 mod 15).  The mirror image bins are then 4(=15-11) and 2(15-13).  Doubling the interval doubles the spacing between the bins.

Side comment:  Fat32's example also shows why I prefer a $1/N$ normalized DFT.  If it were used, the bin magnitudes would be $1/2$ regardless of $N$.

Response to Fat32's comment:

Your approach is still overly cumbersome.  I think it is because you are thinking about the problem backward, and then backsolving.

Step 1. Set Sampling Rate

Step 2. Find N

Step 3. Calculate Interval length (Impicit)

For your first case:

Step 1. $F_s=4000$

Step 2. $N=80$

Step 3. $T_{frame} = 1/50 = .02$ seconds

For your second case:

Step 1. $F_s=375$

Step 2. $N=15$

Step 3. $T_{frame} = 1/25 = .04$ seconds

Clearly your second case frame interval is twice as long as your first case.  Neither your 375 sample rate or the 4001 sample rate will result in an interval that is the .02 second minimum, thus the 28 and 29th bin stipulation is not applicable.

In contrast, I did it in reverse order to yours:

Step 1. Find Interval

Step 2. Set N

Step 3. Calculate Sampling Rate

The minimum frame interval for which both tones have a whole number of cycles is independent of the sampling rate.  Thus, it can be found before any step into the discrete domain needs to be taken.  Finding the minimum value of .02 seconds was done above and I don't think too difficult to understand.  The only intervals that lead to leakage free DFTs are whole integer multiples of the minimum value.  Any non-integer multiplier will lead to non-whole cycle tones.

Once an interval has been chosen, it quite easy to set the number of sample values (N) to any value you want.  It has to be an integer by definition.

The DFT doesn't care about the sampling rate, it only cares about the number of samples.  The sampling rate is calculated from the frame interval duration and the number of samples.

The most important concept is that the bin index (zero based, not MATLAB's stupid one based) corresponds to the number of cycles per frame.  Thus, no matter what N is, as long as N is greater than 29, it will be bin 28 and 29 that have non-zero values for any DFT done on a .02 second interval of this signal.  The value of N will definitely determine where the mirror image bin values are at.  In your second case, since then frame interval is twice as long, the two tones will have frequencies of 56 and 58 cycles per frame respecively.  So those are the bin indexes (or modulus N thereof) that will have non-zero values along with their mirror images.

All possible solutions can then be found by:

Step 1. $T_{frame} = .02 \cdot m $

Step 2. Set an $N$

Step 3. $F_s=N/T_{frame}$

Where $m$ and $N$ are integer values.

Response to Fat32's followup comments:

Comment 1: You are absolutely correct.  Poor wording on my part.  What I meant was "determined" not "constrained".  I have fixed it.

Comment 2: The OP said "The task is to chose the sampling frequency and the base interval of the discrete signal so that there will be no leakage-effect after calculating the DFT."  I think the salient concept behind the question is that you need to have a whole number of cycles of the constituent tones in order not to have leakage in the DFT.  The sampling rate/number of samples is not important in this regard.  I never called your approach wrong, just cumbersome and overly complicated.

Technically, your approach does not work for any $F_s$.  Integer yes, even rational yes, but irrational no.  There are an uncountable infinity worth of values that won't work compared to a mere countable infinity that will.

I think this is poorly worded:  "the resulting N will take more samples than base interval minimum." Perhaps "the resulting N will take a longer interval than the minimum interval."

C3: I did cover the possibility and consequences in my answer of the bin values falling in the DC or Nyquist bins for some choices of N.  Per the criteria of the problem statement, what is important is that there will still not be leakage if this happens.

I am not a MATLAB user, but I am aware of the one based indexing and am also aware that it is the most common criticism of MATLAB.  If it is the tool of choice for you, fine, but arguing what is the best platform quickly devolves into a religious war and is supposed to be off-topic for SE.

SUMMARY:  The problem wasn't stated in either of those ways.  If either additional criteria were given then that would be a discriminant on which solution to choose.  Absent that, I believe my solution demonstrates the underlying concepts clearer and definitely requires fewer calculations.


Measure of harmonicity in a time-series I'm analyzing speech signal for identifying voiced and unvoiced regions. Voiced regions are supposed to have a "pitch", which can be estimated using auto-correlation function (ACF). Basically, one estimates ACF for each frame of speech (say 20ms) and then finds the time lag between peaks in ACF output. If all the significant peaks in ACF output are equidistant, I can say that the signal is very much periodic. If the peaks have random spacing between the, that would indicate aperiodicity.

Based on this, what measure can I use to find out HOW MUCH periodic a signal is? If I denote time lag between successive peaks as lag1,lag2...lagN; the deviation from mean of these values can tell how much periodic the signal is. Any better ideas?

Using the standard deviation of time lags is not a bad idea - the problem is that for very noisy signals such as consonants you won't really get a pattern with peaks. Your suggestion would be more useful in the context of musical instruments sound (for example to measure the inharmonicity of a sound, from violin to piano to bell...)

You can look at the ratio between the value taken at the highest peak in the vocal range (say 80-400 Hz), and the value of the ACF at 0 - this is a common criterion used for voiced/unvoiced classification. Another simple metric that works relatively well with clean recordings is the zero-crossing rate.


How to get the orientation value of a rotated square? Suppose we have a square shape rotated by $d$ degrees counterclockwise (or $-d$ degrees clockwise) in a binary image as the following:

I'm trying to rotate it by 'Orientation' value of regionprops using imrotate function.

But in this case, the major and the minor axis are the same (the ellipse here is a circle).

Hence, the Orientation value equal to $0$ degrees instead of $d$ (or $-d$) degrees.

What is the easiest way to rotate it so that its sides become parallel to the x-axis and y-axis?

(i.e. How to get the value of $d$?)

Disclaimer:  I am not a MATLAB user, so I had to search for this and have no way to test it.

You should be able to find two opposite corner points by using the "Extrema" property from "regionprops".  The tangent of an angle is a slope, so calculate the slope and take the inverse tangent.  Do note that the line may be vertical, i.e. undefined slope, and you have to deal with that case.  If you choose the left-top and right-bottom, you should get a more horizontal line and not have to worry about it.

You also have to keep radians and degrees straight.

As math:

$$ \Delta y = y_{LT} - y_{RB} $$
$$ \Delta x = x_{LT} - x_{RB} $$

$$ m = \frac{\Delta y}{\Delta x} $$

$$ \theta = tan_{deg}^{-1} (m) $$

$$ d = \theta - 45 $$

By finding the leftmost and rightmost $x$'s, there shouldn't be a problem with the slope calculation.  The line between the two extreme points should be a diagonal of the square.  A diagonal should be at 45 degrees, so $d$ is how much you need to rotate the figure to get it there.


I can't wrap my mind around the concept of filtering and the unit circle ? any help? Can someone elaborate this theory where there is a link between the stability of a filter and the zeros or nulls being on the unit circle?

BIBO stability of LTI systems implies that their impulse response is absolutely summable, that is,
\begin{equation}
\sum_{n=-\infty}^{+\infty}|h(n)| < +\infty
\end{equation}

That exact same relationship is a sufficient condition for the Fourier Transform of the impulse response - the so-called Frequency Response - to converge.

Convergence of the frequency response means that it takes finite values, that is,
$$|H(e^{j\omega})|<+\infty \Longrightarrow \sum_{n=-\infty}^{+\infty}|h(n)|<+\infty$$

The frequency response is the evaluation of the Z-Transform on the unit circle: $$H(e^{j\omega})= H(z)\Big|_{z=e^{j\omega}}$$

Moreover, a pole is a point on the complex plane where the Z-Transform's magnitude $|H(z)|$ takes an infinite value: if $z=z_p$ is a pole, then $|H(z_p)| = +\infty$.
Given that the Fourier Transform is a "version" of the Z-Transform on the unit circle, if there is a pole on the unit circle, the frequency response does not converge. If it does not converge, the system is not BIBO-stable.

The last two statements hold because
\begin{equation}
\sum_{n=-\infty}^{+\infty}|h(n)| < +\infty
\end{equation}
does not hold.


Algorithm for testing input sound for multiple expected notes? ** originally I asked this question on stack exchange and was pointed to dsp.stackexchange **

I need to determine when the correct musical notes (from a piano) such as in a chord are being played. While detecting single notes is a no brainer, and full music transcription seems to still be considered a research problem, I'm hoping that since my program knows which combination of notes to expect and only needs to classify the played notes (maximum of six) as correct/incorrect that the task will fall somewhere in between.

Can anyone point me towards an efficient algorithm (I need to implement it in JavaScript) that should be able to achieve this with minimal false negatives (incorrectly telling the user they hit the wrong note)?

My current idea is to test that the signal has better correlation with the target frequency than the nearest semitone above and below for each expected note.

I'd like to point out two things beforehand:

this is a complex question, and
as soon as you have figured out how to do this algorithmically, you'll have to write a lot of signal processing functionality in JS that other languages would already have access to (e.g. python through scipy/numpy, C/C++/Fortran through LAPACK/BLAS/FFTW etc)

so this will become a very complex project, and if you have the chance to influence the expectations for fast results, I'd recommend not promising a working prototype too soon.

Now, for a start, this looks like a linear algebra problem. That's because LA is kind of the sonic screw driver of DSP; give a man a sonic screwdriver, and everything looks like... wait, that metaphor doesn't work.

So, anyways, what you want to do is find a vector representation of your signal under test and map it to coordinates in a space spanned by the possible musical notes. Mathematically, you want to

Find a set of base vectors representing your notes,
Find a suitable base in which you can represent both these vectors and your signal
Project the signal vector onto your note vector; the magnitude of the coefficients gives you the "amount" of the respective notes being represented in the signal.

As an example, the discrete Fourier transform takes signals that are represented in a time sampling base (every base vector is $(0,\ldots,0,1,0,\ldots,0)$, with the 1 representing the the sample time) and transforms it to a base system of orthogonal oscillations (The nth base vector of that base is $( e^{j2\pi nt})_{t\in (0,\ldots, N)}$). That is such a base transform.

You could simply record fixed lengths of your tones, normalize them in power (important!) and do a dot product with each of these, and find the ones that are above a certain threshold. Do that again with the piece of your signal starting one sample later, and find the notes that have been present with enough energy for enough time. What you end up with is a FIR filter bank.

Now, with a proper language and libraries, that would be a few lines of code, and one could do everything in frequency domain rather than time domain, reducing the effort of the convolution above to a multiplication and two FFTs, but since there is (I don't think there can be, due to architectural restrictions) no FFT in JavaScript, this is probably not easy to optimize. It might be that for reliable results, your algorithm might become too slow to process the audio signal in real time.


Why does stability only consider, whether the system is bounded? That doesn't mean the response is smooth? Why does (BIBO) stability only consider, whether the system is bounded? That doesn't mean the response is smooth?

That is, one could have system whose values are bounded. Still it could exhibit "chaotic" behaviour through its parameters.

So what is this kind of stability/instability concept really helpful for?

BIBO stability applies to linear systems with inputs and outputs.  See

https://en.wikipedia.org/wiki/Lyapunov_stability

There is a section on linear systems with inputs.  The article is clear and restating it here is unnecessary.


machine learning on multiple signal at once I am very new signal processing and I am asking this with the expectation that someone gives me a concrete guideline on how to approach my problem.

The problem I need to solve is that I have 2 signals (shown in green and yellow) which look like this as below in original scales

and log scale

Based on the interaction (e.g. cutting each other) of these 2 signals, few features are calculated like i.e. range where both signals are low, the range where signals are up and only a threshold difference in amplitude value but now I need some kind of Machine learning may be RNN or LSTM or 1D CNN. I have a set of pre-calculated values and corresponding signals. Now I am wondering how to apply machine learning techniques while looking at both these signals together.

You say you're "a beginner" so I will suggest using Keras (https://keras.io/). It is written in Python and runs on top of Tensorflow which is a neural network toolbox written by Google. Keras allows you to build a network layer-by-layer in an easy way. For example, here is an example of building a 1D CNN (https://keras.io/examples/imdb_cnn/), please note that the example contains many things that you do not need or are optional. My advice, start with a simple network and know how to identify its problems and how to fix them.

A while back, I worked on a similar problem involving inputting two time domain signals and classifying it (input: [# samples]-by-2, output: probabilities corresponding to each class). Now you will input the two time domain signals and you need to think carefully about how your output is formatted because sometimes you might have one activity, other times you might have five activities, or sometimes you might have no activity.

For example, say you input two signals ($x_1[n], x_2[n]$) of length 1024 and say you had activity #1 during samples 0-255, activity #2 from 256-511, and all other time is activity #0. Then you could format your output to be a signal of length 1024 (equal to length of input signals), to be:
$y[n]=
 \begin{cases}
      1, & 0 \leq n \leq 255 \\
      2, & 256 \leq n \leq 511 \\
      0, & n \geq 512
   \end{cases}
$. I bet this requires a bit of preprocessing of your dataset but this is not uncommon at all when doing these sorts of things.


Inverse Discrete Time Fourier Transform of $1$ $\textrm{DTFT}(\delta[n]) =1$,
but $\textrm{IDTFT(1)} = \frac{\sin(\pi n)}{\pi n}$.
Why it is not equal to the unit impulse $\delta[n]$?

The IDTFT of $X(e^{j\omega})=1$ is indeed

$$x[n]=\frac{\sin(n\pi)}{n\pi}\tag{1}$$

Now, what happens for indices $n\neq 0$? As it turns out, you can safely rewrite $(1)$ as

$$x[n]=\delta[n]\tag{2}$$

where $\delta[n]$ is the discrete-time unit impulse. (HINT: think about where the zeros of $\sin(x)$ are).


Why not use the same "standard" exponentials for both continuous and discrete time In continuous time the standard exponential signal is usually defined as
$$
e^{st}, \quad\text{with}\quad s = \sigma+j \omega
$$
In discrete time the standard exponential signal is usually defined as
$$
z^n, \quad\text{with}\quad z=re^{i \theta}
$$
Why not define the discrete time signal as $ e^{zn} $? It seems that $ e^{zn} $ is also an eigenfunction for discrete LTI systems.

This has to do with the way the Laplace transform and the $\mathcal{Z}$-transform are defined:
$$\mathcal{L}\big\{x(t)\big\}=\int_{-\infty}^{\infty}x(t)e^{-st}dt\tag{1}$$
$$\mathcal{Z}\big\{x[n]\big\}=\sum_{n=-\infty}^{\infty}x[n]z^{-n}\tag{2}$$
Of course, the question remains why they are defined that way. The answer lies in the fact that in continuous time we usually have to deal with differentiation and integration, which in the Laplace transform domain correspond to multiplication with $s$ and $1/s$, respectively. Whereas in discrete time, we usually have to deal with difference equations formulated in terms of time delays and time advances, which in the $\mathcal{Z}$-transform domain correspond to multiplications by $z^{-1}$ and $z$, respectively.
Consequently, the definitions $(1)$ and $(2)$ lead to transfer functions of linear time-invariant (LTI) systems that are ratios of polynomials in $s$ for continuous time, and polynomials in $z$ for discrete time. In this way, the choices of $s$ and $z$ in $(1)$ and $(2)$ and the corresponding definitions of exponentials as $e^{st}$ and $z^n$ appear natural.


tx and rx model supporting all data rates i am trying to build complete tx and rx path which can support varies data rates. for example if my data rate have to support from 10 kbps to 1000 kbps.
Can i keep the baseband processing fixed for all the data rates same. what i mean is using a CIC filter in tx and rx which always results in processing
4 samples/symbol. or should i use a resampler in the timing recovery which will result in the desired rate?

In general terms, it is usually advantageous to down-sample to the lowest sampling rate wherever possible- all things being equal it would be better to reduce to a lower sampling rate earlier in the receiver processing - for purposes of reducing resources required, power dissipation, etc.
The OP is using both a CIC and timing recovery block where the opportunity to change the symbol rate can be done in either location (for example to 4 samples per symbol in any case). Here the preferred approach would likely be in the CIC, to the extent that all subsequent blocks can run more efficiently, with the same performance metrics, at the lower sampling rate. If any such block would require or be improved with a higher sampling rate, then resampling further downstream would be preferred. This would be clear in the detailed implementation and specific to that. Also note that a CIC introduces a passband droop that may need to be compensated for. This can be easily implemented with a 3 tap FIR filter as I have demonstrated at this post.
how to make CIC compensation filter


How to add white gaussian noise to an image in Matlab? I want to add white gaussian noise to an image of 10 dB in Matlab.

I tried to use Matlab function imnoise but I couldn't figure out what values for mean and variance should I choose to add noise of 10 dB.

How do I do this ?

Additive noise is usually mean-free, i.e. the mean should be 0. A noise level of 10 dB = 1 B ("bel") usually means that the variance of the noise is by a factor 10¹ = 10 smaller than the variance of the image. As the documentation states, the variance parameter to imshow is interpreted under the assumption that the image data are doubles in the range [0, 1].

Here's an example:

Load sample grayscale image included with Matlab:

[I, map] = imread('eight.tif');

That image is of data type uint8, i.e. integers in the range 0 to 255. Let's convert it to a double in [0, 1]:

I = double(I) / 255;

Compute the image variance

v = var(I(:));

and tell imnoise to add Gaussian white noise with mean 0 and a tenth of that variance:

J = imnoise(I, 'gaussian', 0, v / 10);

This is the result:

Because imnoise clips the data to the range [0, 1] after adding noise, the actual noise variance is smaller,

var(J(:) - I(:)) / v

gives something like 0.093.


Logarithm in state space equations I want to linearize a system to this form
$$\begin{bmatrix}
   \Delta\dot{x}_1\\
  \Delta\dot{x}_2\\
  \Delta\dot{x}_3
 \end{bmatrix} = A\begin{bmatrix}\Delta x_1\\ \Delta x_2\\ \Delta x_3\end{bmatrix}+B\begin{bmatrix} \Delta u\end{bmatrix}$$
    $$\begin{bmatrix}
   \Delta y
 \end{bmatrix} = C\begin{bmatrix}\Delta x_1\\ \Delta x_2 \\ \Delta x_3\end{bmatrix}+D\begin{bmatrix} \Delta u\end{bmatrix}$$

The A matrix

$$A = \begin{bmatrix} \frac{\partial \dot x_1}{\partial x_1} & \frac{\partial \dot x_1}{\partial x_2} & \frac{\partial \dot x_1}{\partial x_3}\\[6pt] \frac{\partial\dot x_2}{\partial x_1} & \frac{\partial \dot x_2}{\partial x_2} & \frac{\partial \dot x_2}{\partial x_3}\\[6pt] \frac{\partial\dot x_3}{\partial x_1} & \frac{\partial \dot x_3}{\partial x_2} & \frac{\partial \dot x_3}{\partial x_3}\end{bmatrix}_{|P}$$

$$=\begin{bmatrix} M_1L & \frac{d}{dt} \left[\frac{R_2}{M^2}\ln(x_2(t)+1)\right] & R_1\\[6pt] L_1L_2   &   \frac{d}{dt} \left[\frac{L_1}{M^2}\ln(x_2(t)+1)\right]& R_2\\[6pt] 0 &\frac1C & 0\end{bmatrix}_{|P}$$

has in its second comulmn logarithm. Pay attention only there. ($M$, $L$'s and $R$'s are some constants) If I derivate it, I get

$$=\begin{bmatrix} M_1L & \frac{R_2 \frac{dx_2(t)}{dt}}{M^2(x_2(t)+1)} & R_1\\[6pt] L_1L_2 & \frac{L_1 \frac{dx_2(t)}{dt}}{M^2(x_2(t)+1)}& R_2\\[6pt] 0 &\frac1C & 0\end{bmatrix}_{|P}$$

Am I right?

And I'm not sure what to do with this. Since my initial condition $P$ are in equilibrium and  $\frac{dx_2(t)}{dt} = 0$, should I put there the $0$? The result would look like

$$=\begin{bmatrix} M_1L & 0 & R_1\\[6pt] L_1L_2 & 0 & R_2\\[6pt] 0 &\frac1C & 0\end{bmatrix}_{|P}$$

But few weeks ago I have been solving similar situation and I didn't put $\frac{dx_2(t)}{dt}$ in the matrix at all. Lecturer told me it was right. According to that, this matrix should look like

$$=\begin{bmatrix} M_1L & \frac{R_2}{M^2(x_2(t)+1)} & R_1\\[6pt] L_1L_2 & \frac{L_1 }{M^2(x_2(t)+1)}& R_2\\[6pt] 0 &\frac1C & 0\end{bmatrix}_{|P}$$

But I'm confused. The mode $x_i(t)$ isn't consider as a function of time?

How to solve this?

It's $\frac{d\dot x_2}{dx_2}$, not $\frac{d\dot x_2}{dt}$, so that's why the right result truly is $\frac{R_2}{M^2(x_2(t)+1)}$.


Group delay of the FIR filter For an FIR filter, with symmetrical tap values $h[N-1-n]=h[n]$,

why is the group delay $\frac{N-1}{2} T$ (where $N$ is the number of taps of the FIR filter and $T$ is the sampling time)?
Why is linear phase so important for the filter response?

To be precise the group delay of a linear phase FIR filter is $(N-1)/2$ samples, where $N$ is the filter length (i.e. the number of taps). The group delay is constant for all frequencies, because the filter has a linear phase, i.e. its impulse response is symmetrical (or asymmetric). A linear phase means that all frequency components of the input signal experience the same delay, i.e. there are no phase distortions. So for a frequency selective filter (e.g., a low pass filter), if the input signal is in the passband of the filter, the output signal is approximately equal to the input signal delayed by the group delay of the filter. Note that in general FIR filters do not have a linear phase response. In this case, the group delay is a function of frequency.

EDIT:

Consider an odd length FIR filter with even symmetry ("type I filter"). In this case the group delay is an integer number of samples $K=(N-1)/2$. The filter has a linear phase response (because the group delay is the negative derivative of the phase) and the filter's frequency response can be written as

$$H(\omega)=A(\omega)e^{-jK\omega}$$

where $A(\omega)$ is the real-valued amplitude function, and $\phi(\omega)=-K\omega$ is the (linear) phase. If the input signal's Fourier transform is $X(\omega)$, then the Fourier transform of the output signal is given by

$$Y(\omega)=H(\omega)X(\omega)=A(\omega)X(\omega)e^{-jK\omega}\tag{1}$$

Now assume that the frequency content of the input signal is in the passband of the filter, i.e. in the frequency region where $A(\omega)\approx 1$ holds. Then from (1) the output signal's spectrum can be written as

$$Y(\omega)\approx X(\omega)e^{-jK\omega}\tag{2}$$

and, from the shift property of the Fourier transform, the output signal is (approximately) given by

$$y(n)\approx x(n-K)$$

i.e. $y(n)$ is a delayed version of the input signal, where the delay is given by the group delay $K$.


What is this theorem in this formula? I just read a paper yesterday and there is confusion:

I want to know what is reconstruction fidelity term and the theorem of changing the formula from upper to lower one. Is this theorm from functional analysis?

Could anyone help me ?

EDIT: The link of the paper.

i have no idea what the "reconstruction fidelity term" is or what it's about.

Hermitian symmetry is a term usually applied to some form the Fourier Transform of a signal that is purely real.

for continuous-time, continuous-frequency Fourier Transform:

$$ X(f) = \int\limits_{-\infty}^{\infty} x(t) \, e^{-j 2 \pi f t} \ dt $$

if $x(t)$ is purely real (that is $\Im\{x(t)\}=0$ for all real $t$), then we know that there is this symmetry about $f=0$:

$$ X(-f) = X(f)^* = \operatorname{conj}\{X(f)\} \qquad \forall f \in \mathbb{R} $$

or

$$\begin{align}
\Re\{X(-f)\} &= \Re\{X(f)\} \\
\Im\{X(-f)\} &= -\Im\{X(f)\} \\
|X(-f)| &= |X(f)| \\
\arg\{X(-f)\} &= -\arg\{X(f)\} \\
\end{align}$$

Similarly, for the Discrete Fourier Transform:

$$ X[k] = \sum\limits_{n=0}^{N-1} x[n] \, e^{-j 2 \pi \frac{nk}{N}} $$

if $x[n]$ is purely real (that is $\Im\{x[n]\}=0$ for integer $n \in [0..N-1]$), then we know that there is this symmetry about $k=\tfrac{N}{2}$:

$$ X[N-k] = X[k]^* = \operatorname{conj}\{X[k]\} \qquad 1 \le k \le N-1 $$

or

$$\begin{align}
\Re\{X[N-k]\} &= \Re\{X[k]\} \\
\Im\{X[N-k]\} &= -\Im\{X[k]\} \\
|X[N-k]| &= |X[k]| \\
\arg\{X[N-k]\} &= -\arg\{X[k]\} \\
\end{align}$$

Now in my religion regarding the DFT, i insist that

$$\begin{align}
x[n+N] &= x[n] \qquad & \forall n \in \mathbb{Z} \\
X[k+N] &= X[k] \qquad & \forall k \in \mathbb{Z} \\
\end{align}$$

this is always the case (the DFT maps a discrete and periodic sequence of period $N$ in one domain to another discrete and periodic sequence of period $N$ in the reciprocal domain).  then the Hermitian symmetry takes a simpler form:

if $x[n]$ is purely real (that is $\Im\{x[n]\}=0$ for all integer $n$), then we know that there is this symmetry about $k=0$:

$$ X[-k] = X[k]^* = \operatorname{conj}\{X[k]\} \qquad \forall k \in \mathbb{Z} $$

or

$$\begin{align}
\Re\{X[-k]\} &= \Re\{X[k]\} \\
\Im\{X[-k]\} &= -\Im\{X[k]\} \\
|X[-k]| &= |X[k]| \\
\arg\{X[-k]\} &= -\arg\{X[k]\} \\
\end{align}$$

likewise you can identify Hermitian symmetry for the DTFT and relate it to the real-ness of the input to the DTFT. (do i have to do that one, too?)


How to Calculate the Color Temperature / Tint of the Colors in an Image? How can I calculate the (black body equivalent) color temperature of the light source illuminating a given image? Below are screenshots of Adobe Lightroom manipulating color temperature and a shift in the RGB Histograms. Given RGB components of the image, how do I go about calculating it? I should expect a single value - The black body equivalent temperature of the illumination source, right?

This paper (pdf download) gives the following formulae for calculating Correlated Color Temperature (CCT). They do not explicitly say (or I missed it), but their example leads me to infer that they are assuming RGB values in the range of 0-255.

1. Convert the RGB values to CIE tristimulus values (XYZ) as follows:

$$X = (-0.14282)(R) + (1.54924)(G) + (-0.95641)(B)$$
$$Y = (-0.32466)(R) + (1.57837)(G) + (-0.73191)(B) = Illuminance$$
$$Z = (-0.68202)(R) + (0.77073)(G) + (0.56332)(B)$$

2. Calculate the normalized chromaticity values:

$$x = X/(X+Y+Z)$$
$$y = Y/(X+Y+Z)$$

3. Compute the CCT value from:

$$CCT = 449n^3 + 3525n^2 + 6823.3n + 5520.33$$

$$\text{where }n = (x − 0.3320) / (0.1858 − y)$$

Which can be combined to form the following equation:

$$CCT = 449n^3 + 3525n^2 + 6823.3n + 5520.33$$
$$\text{where }n = ((0.23881)R+(0.25499)G+(-0.58291)B) / ((0.11109)R+(-0.85406)G+(0.52289)B)$$

I am not sure about applying this to an image, but if you just want a single, generalizing number, than you could perhaps use some sort of averaging? Either find an acceptable "average" RGB value for the image (eg. the centroid) and use that to calculate a temperature or (a much more computationally expensive option) calculate the temperature for each pixel in the image and take the average of those results.

Also, bear in mind that CCT is only an approximate metric for most colors, since only a single curve in the color space actually represents color which can be obtained from a real world black body radiator. Thus for all other colors, the calculated color temperature is simply an approximation of the black body temperature it most closely represents. Thus, for some colors (especially greens) it can actually be a somewhat meaningless value, at least in a physical sense. This is illustrated well in the following image (from the wikipedia article on color temperature).

The black line in the image represents the Planckian Locus of colors which could actually be produced by block body radiation. The smaller crossing lines represent the isotherms of the CCT approximation nearbye.

Also, since your question specifically references Adobe Lightroom, I found this while searching around:

  The sliders [in Adobe Lightroom] adjust not the black body temp of the light, but the compensation applied to the image to compensate for the black body temp of the light.  This goes the other way round.

So bear in mind that the color temperature you see on the Lightroom slider will not be the same as those calculated from the above formulae.


How to interpret spectrogram correctly? I have a pressure transducer signal from a reacting mixture available here, which has a sampling frequency of 50kHZ. I have also obtained the spectrogram of the signal with NFFT=256, noverlap=4 as shown below:

I expected to see one high power frequency after pressure rise, instead of multiple frequencies (after t=0.03s). Why are multiple frequencies detected? Am I using the correct NFFT and noverlap parameters? How are NFFT and noverlap selected for this signal?

The problem is not the spectrogram parameters, these are correct since they only depend on what resolution you want in time and frequency domain. Also, the spectrogram interpretation is correct, there are multiple frequency peaks. The problem may be:

  I expected to see one high power frequency after pressure rise, instead of multiple frequencies

Why? If you are not in free-field conditions (and you mention a cylinder) there will be multiple wave reflections, resonances, standing waves and 3D pressure patterns. A sudden pressure step will distribute its power over the spectrum, but frequency peaks will appear because you are essentially "tuning" your cylinder.

Recall that pressure waves travel at sound speed. Considering for example the habitual 340 m/s; if your cylinder measures 1m, the wave front will travel back and forth each 0.0029s, adding multiple frequency components to the measurement.

As for standing wave patterns, we can consider the case of a closed cylinder, where we can compute the frequency of the fundamental frequency as $f_1 = \frac{a}{4L}$ with subsequent harmonics at $nf_1$ where $n$ are odd integers (if your cylinder is indeed closed, does it have a length of ~20cm?).

Regarding the problem of measuring higher frequency waves in ducts, consider that according to Eriksson, it may be assumed that waves in a circular duct propagate in a planar fashion (i.e., pressure is constant in every cross-sectional area of the duct) only below a cut-off frequency of:
$$
f = 1.84 \frac{a}{\pi D}\sqrt{1-M^2}
$$

Where a is the sound speed, D the duct diameter and M the Mach number of the mean flow. Above this frequency, the pressure distribution in a given section of the duct will not be uniform, it will have nodes. See the following animation for an idea on how pressure varies in a given planar section, noting that only the Z(0,0) is a typical planar wave.

Beyond the planar wave frequency it is difficult to correctly measure the pressure spectrum with only one transducer. In any case, I believe that the signal processing that you are performing is correct; the only problem is that you are witnessing the acoustic effects of your cylinder.


FIRFilter with low cut-off compared to sampling frequency Sorry if my question is some FIRFiltering basics, but I couldn't find a specific answer after some search.

I am trying to design a decimating low pass filter in C (actually in OpenCL with GPU) with the FIR approach by doing convolution.
I have a software that generates the coefficients based on the cut-off frequency and the sampling frequency.

I realized that my filter generally works if the cut-off frequency is not too far away from sampling frequency, for example fsample = 1e7 and fcutoff = 1e4 would work. (maybe the attenuation is not enough). But if I increase the fsample to 1e8, the attenuation will be negligible.

Reading from this link:
How many taps does an FIR filter need?
I understand that I will need a number of tap of order 1e5 or 1e6 to have sufficient attenuation (~10 - 20 dB) for my application.

So my question is, in general, how should I design my filter to meet my requirement, if I really have to cut-off something massively oversampled?
If not FIR, could you point me any other better approach, for example, can IIR work better?

Maybe I can follow the approach of multi-stage decimating FIR filter in this page? Great if you provide me with some references regarding how the stages should be determined.
FIR-Decimation and Low-pass filter (taps vs number of input points vs number of decimation stages)

Thanks.

If the desired cutoff frequency is very low relative to the sample rate, then yes, you would need a very high-order FIR filter to implement it. However, in these sort of cases, it's usually a better idea to use a multistage approach.

For example, say you wanted to decimate a signal by a factor of 10000. If you implement this with a decimating lowpass filter, it would need to have a cutoff frequency less than $\omega = \frac{\pi}{10000}$. This would require a filter with tens to hundreds of thousands of taps, depending on how much attenuation you require in the stopband. This is a lot of processing, especially if you're trying to do real-time processing. What's more, with really long filters like this, you can run into numerical precision issues; the filter that you get might not even give you the frequency response that you designed it for!

Instead, consider an approach with two stages: two identical filters that each decimate by a factor of 100. In this case, your cutoff frequency must be less than $\omega = \frac{\pi}{100}$. And you'll only need hundreds to thousands of taps in order to implement each filter. The computational cost of two few-thousand-tap filters is much less than the cost of one filter with hundreds of thousands of taps. This only becomes more apparent when you notice that the first filter decimates by a factor of 100, making the cost of the second one negligible compared to the first stage.

In short, when large rate changes like this are desired, multi-stage approaches are almost always what you want. The same logic follows for interpolation as well. Lyons has an easy-to-understand discussion of this topic in his introductory text Understanding Digital Signal Processing.


white noise in images For digital images, noise is assumed to be additive gaussian white noise. I remember noise in images is considered as high frequency. However, from the link
https://en.wikipedia.org/wiki/White_noise, it says

  a random signal is considered "white noise" if it is observed to have
  a flat spectrum over the range of frequencies that is relevant to the
  context.

Does it mean frequencies of noise will extend from low frequencies to high frequencies? Why do we try to only remove high frequency component? Is human visual system  not sensitive to low frequency noise?

From the same link, https://en.wikipedia.org/wiki/White_noise

  In digital image processing, the pixels of a white noise image are
  typically arranged in a rectangular grid, and are assumed to be
  independent random variables with uniform probability distribution
  over some interval.

It reads that the pixels of a white noise image are assumed to be independent random variables with uniform probability distribution over some interval.
Does it mean that one pixel and its neighbors are i.i.d?

Yes, the noise is wide-band.

A challenge in images resides in lowering the high-frequency noise components, while preserving the local high-frequency features of the contours and the textures, where you have fewer pixels to estimate a consistent denoised value.

On low-frequency parts, the data is generally more regular or smooth, and getting rid of the noise is often performed more easily by averaging over neighboring pixels.


How can I measure the harmonicity of an audio spectrum? I would like to estimate the "harmonicity" (harmonics-to-noise ratio ?) of an audio signal from its spectrum. What kind of algorithm can I use?

Please take a look at this previous answer.  Using the semantics and mathematical definitions from that answer...

Let $x[n]$ be the input signal from which you want to measure harmonicity or periodicity.  Assume $x[n]$ has no DC component (like it is the output of a DC blocking filter).  And $n_0$ is the sample index of the neighborhood of signal around where you want the harmonicity to be measured.  $N$ is a large integer defining how wide your rectangular window is.

Average Squared Difference Function, ASDF:

$$ Q_x[k, n_0] \triangleq \frac{1}{N} \sum\limits_{n=0}^{N-1} \left(x[n+n_0-\left\lfloor \tfrac{N+k}{2}\right\rfloor] \ - \ x[n+n_0-\left\lfloor \tfrac{N+k}{2}\right\rfloor + k] \right)^2 $$

$\left\lfloor \cdot \right\rfloor$ is the floor() function and, if $k$ is even then $ \left\lfloor \frac{k}{2}\right\rfloor = \left\lfloor \frac{k+1}{2}\right\rfloor = \frac{k}{2} $.

From that define an alternative form of Autocorrelation:

$$ R_x[k,n_0] = R_x[0,n_0] - \frac12 Q_x[k, n_0]  $$

where

$$ R_x[0, n_0] \triangleq \frac{1}{N} \sum\limits_{n=0}^{N-1} (x[n+n_0-\left\lfloor \tfrac{N}{2}\right\rfloor])^2 $$

Since $Q_x[0, n_0] = 0$ and $Q_x[k, n_0] \ge 0$ for all lags $k$, that means that $ R_x[k, n_0] \le R_x[0, n_0] $ for all lags $k$.

You will find a period $P$ such that $P > 0$ and $R_x[P, n_0] > R_x[k, n_0]$ for all lags $k>k_0$, other than the highest peak at $k=0$.

$k_0$ is the smallest value of $k>0$ such that $ R_x[k, n_0]<0 $.

Once that peak location $P$ is found, the harmonicity is

$$ r_x(n_0) \triangleq \frac{R_x[P, n_0]}{R_x[0, n_0]} $$

If $r_x(n_0) = 1$, then $x[n]$ is perfectly periodic (or "harmonic") in the vicinity of sample index $n_0$. If $r_x(n_0) = 0$, it's noise, no periodicity to be found.