Frequency Response Notation Why do people sometimes write $H(e^{j\omega})$, some others use $H(\omega)$ and even others use $H(f)$ with $\omega=2\pi f$ to describe the frequency response of a filter or a spectrum? Is there any "real" difference in meaning or is it just convention?

In nearly every field of mathematics or science you will find multiple notations that mean the same thing, sometimes from slightly different angles. This is one of those cases where there are multiple ways of notating the same concept.

For discrete-time systems, there really isn't any difference in meaning between the definitions that you gave, especially between the first two ($H(e^{j\omega})$ and $H(\omega)$). The first is just a bit more explicit that the function $H$ refers to the (complex-valued) amplitude of the system's output when its input is a complex exponential with angular frequency $\omega$. This is a bit verbose, so often you'll see it trimmed to the latter representation. $H(f)$ is not typically used for discrete-time systems, as $f$ refers to a continuous-time frequency in Hertz; without a specified sampling frequency, there isn't a well-defined mapping between $f$ and $\omega$, the more natural frequency concept for discrete-time systems.

$H(f)$ is typically only specified for continuous-time systems, where it competes with the alternate notation $H(\omega$). In any mathematical expression where you use one of them, you can massage the expression to use the other, sometimes involving factors of $2\pi$ that need to be inserted in various places to make the math work out. As a general rule, in my experience, the $\omega$ form is used more commonly in control systems contexts, whereas communications engineers seem to prefer the $f$ form. However, as I said, there's really no functional difference between the two.


Why does downsampling increase the phase angle of poles in autoregressive models? When reading about autoregressive models (AR), people often "downsample to increase the phase angle of the poles of the corresponding transfer function".

How exactly does downsampling affect the phase angle of the poles?

Decimating a signal $x[n]$ by an integer factor $M$ results in the following $\mathcal{Z}$-transform:

$$Y(z)=\frac{1}{M}\sum_{k=0}^{M-1}X\left(e^{-j2\pi k/M}z^{1/M}\right)\tag{1}$$

If you assume that $x[n]$ was filtered by an anti-aliasing filter before decimation, the transfer function of the decimated filter can be approximated by the term for $k=0$ in (1):

$$Y(z)\approx\frac{1}{M}X(z^{1/M})\tag{2}$$

So if $z_{\infty}$ is a pole of $X(z)$ then $z_{\infty}^{M}$ is a pole of $Y(z)$. Writing the pole in polar coordinates $z_{\infty}=re^{j\phi}$ gives

$$z_{\infty}^M=r^Me^{jM\phi}$$

So the angle of the pole is increased by a factor $M$.


Relationship between DFT input sequence and magnitude Assuming there is a sequence that could look like this:

$$ x[n] = \{1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0,1, 1, 0, 0, 1, 1, 0, 0,1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0\} $$

Using this sequence, I want to perform a DFT manually, the connection between $x[n]$ and $|X[k]|$ is important to me.

Now there is only one problem: my sequence contains 32 values, so the DFT by hand would be very complicated. Is there a faster way to see the connection/relationship between the original signal (sequence) $x[n]$ and the magnitude $|X[k]|$.

One can recognize a period in the sequence, perhaps that could be exploited, the question being, how?

To make it short: Assuming I have a sequence (e.g. the above given), how can I infer the magnitude $|X[k]|$ from it?

I hope my question is understandable so far. Thanks for coming answers!

In general, one can't infer the magnitude of the DFT coefficients from the values of the input sequence if by inference one means a process that is much more amenable to manual computing than just evaluating the DFT directly. In a sense, you are asking for "half" (or more) of the information that the DFT provides without doing much work for it; TANSTAAFL. What can be done easily is to get a (weak) upper bound (in some sense, a trivial upper bound) on $|X[k]|$. We have that
\begin{align}
|X[k]| &= \left|\sum_{n=0}^{N-1} x[n]\exp(-j2\pi kn/N)\right|\\
&\leq \sum_{n=0}^{N-1} \left|x[n]\exp(-j2\pi kn/N)\right|\\
&\leq \sum_{n=0}^{N-1} \left|x[n]\right|\\
&\leq N \max_n |x[n]|.
\end{align}
In typical cases of DFT computations, this is a very weak bound. In the
OP's special case however, one can use the fact that $16$ of the $32$ $x[n]$ are $0$ (and the rest have value $1$) to tweak the bound to
$$|X[k]| \leq 16.$$
Since jomegaaA's calculations show that $X[0]$ does indeed have value $16$ in this instance, no tighter bound can be found in this particular instance.


MATLAB Implementation of Karplus Strong algorithm with filter function? I want to implement following function: y = ksalgrithm(x, alpha, M, Nout) where

x is the input vector with length Nout
alpha=0.99
Nout is the length of the output vector.

The algorithm should do the convolution of $y[n]=(h_1\star x)[n]$, where $h_1[n]$ is the inverse of the $z$-transformation of $H_1(z)$. The input signal should have values $l = T_0 \cdot f$ and filled up with zeros until the vector has the length of Nout.

I want to use MATLAB's filter.

Nout = length(x);
M = T_0*fs;
alpha = 0.99;
function [ y ] = ksalgorithm( x, alpha, M, Nout )
    syms z;
    H = filter(1,1-alpha*z^(-M),x);
    H_1 = 1/(1-alpha*z^(-M));
    h_1 =  iztrans(H_1);

    y = zeros(x);
    ...

I get always this error?

Undefined function 'filter' for input arguments of type 'sym'.
Error in ksalgorithm (line 23)
    H = filter(1,1-alpha*z^(-M),x);

The filter(b,a,x) function of MATLAB does not operate on symbolic variables. Your line of code

syms z

creates a symbolic variable z, in MATLAB environment, and later
you feed this variable into the filter() function through its second argument in this line:

H = filter(1,1-alpha*z^(-M),x);

where the filter(b,a,x) function, in this mode is expecting a numerical coefficient vector $a[k]$ associated with the DT system that's defined by the coefficient vectors $a[k]$ and $b[k]$.

In order to get the filter() function working, replace the symbolic expression

1-alpha*z^(-M)

with the numerical one, given integer M and real alpha

a = [1 zeros(1,M-1) -alpha]


Plotting curved peaks on PSD graph I am a PhD student with little engineering experience, currently attempting to recreate the signal processing techniques used in this paper: Monitoring Vital Signs and Postures During Sleep Using WiFi Signals. As a TL;DR: changes in channel state information (amplitude, in this instance) can be attributed to movements in otherwise still environments, most notably breathing and heart rate.

In recreating the signal processing performed in this paper, I have found it challenging to follow how the authors produced results as seen in Figure 8a and 8b, which should show isolated breathing frequencies.

At this phase, the signal (across 30 subcarriers) has been filtered using:

Hampel filter (for significant outlier removal)
Running mean filter (for smoothing)
Butterworth bandpass filter (to isolate potential breathing frequencies)

Once this has been done, a heatmap of amplitude across the subcarrier spectrum can be produced which should show clear banding, as seen in Figure 8a.
The Power Spectral Density for each subcarrier can be plotted, with the peak of the mean indicating the strongest observed frequency (Figure 8b). However, the lines plotted in this figure are curved, and I have been unable to find another example of curved PSD. It seems as though there is less data than implied in this graph, as this figure represents 5s of data (which was captured at 20Hz), which means the PSD could only be plotted for a 100n FFT. In my understanding, this should mean the bins are 0.2Hz wide? It appears as though they are able to use these curves across the peaks to infer additional information.

Using MATLAB's spectrogram, I have been able to plot the observed PSD in a similar signal (captured at 100Hz), and the peaks on this graph are not curved. My understanding of FFTs is that the peaks for each bin can be observed, but an asymmetrical curve as seen in Figure 8b cannot be drawn as we have no further information other than the bin.

Is my understanding of FFTs or PSD incorrect, or have I missed an obvious step in producing the results as seen?

Different methods exist for estimating PSDs. They can be broadly classified as parametric or non-parametric methods. Indeed you should learn which method they used to obtain the associated PSD.

An FFT based PSD estimation is a non-parametric method such as Periodogram and its variants. This method will produce random looking non-smooth output result which can be smoothed, but this has limits and has side effects too.

On the other hand, parametric estimators can produce very smooth and high resolution looking PSD estimates. There are many parametric PSD estimators. The authors probably used a model based parametric estimator to obtain the plotted PSD, that's why they are smooth looking and that's what I understand from your word curved ?


Do I need to bandpass filter the signal before applying the Zoom FFT algorithm? Let's say I have an acoustic signal that I want to make a spectrogram of filled with broadband noise and tonals, sampled at oh, $44100\textrm{ kHz}$, and I want to zoom in on the band between $2600\textrm{ Hz}$ and $3000\textrm{ Hz}$ (or any band of frequencies, really).

The zoom FFT algorithm says to accomplish this, one simply has to frequency translate, then lowpass filter, then decimate and FFT the resulting lower frequency signal.

However, most of the diagrams I see explaining zoom FFT only show discreet bands with no noise in between, which begs the question:

Would I need a bandpass filter to prevent aliasing during the frequency translation step of the Zoom-FFT algorithm if I have a noisy signal containing energy at all frequencies?

I'm thinking yes, since there's a chance that during the frequency translation some unwanted noise would get shifted into my new band of interest due to aliasing, but I'm not sure.

The processing chain I envision is as follows:
$$\boxed{\text{BPF}}{\rightarrow}\boxed{\text{Frequency Xlate}}{\rightarrow}\boxed{\text{Bandlimited downsample}}{\rightarrow}\boxed{\text{FFT}}{\rightarrow}\boxed{\text{Display processing}}$$

Would this be correct?

Your BPF is superfluous. Frequency translation (if very simply done by $\cdot e^{j2\pi \Delta f n}$) is a circular shift of the whole Nyquist band.

When you calculate the Fourier of that, you'll notice there's no chance of aliasing at all.

The only step that leads to aliasing is the downsampling, and that's why ever decimator comes with an appropriate filter - in your case, the LPF you already have. You usually implement the filtering and down sampling as one entity - no use in calculating filter output that you'll throw away in the next step.

So:

$\cdot e^{j2\pi \Delta f n} \rightarrow \text {decimating LPF}\rightarrow \text {DFT}$


What is the best way(In terms of speed) to implement a direct form type 2 IIR Filter in VHDL? I am trying to implement a 4th order direct form type II IIR filter in VHDL(Specifically, in xilinx FPGA) and when i implement it i get correct output yet very low maximum clock rate(Around 80 Mhz).
So i wanted to know what is the best way to implement IIR filters in terms of speed?Or at least where can i read more about IIR filter implementations in VHDL?

I faced the same problem in the past, the problem is that it's not possible to pipeline the IIR feedback loop without changing the filter equation.
There's one technique you can use, it's called "scattered-look ahead IIR filter"
The gist of it is to change $H(z) = \frac{b_0 + b_1z^{-1} + b_2z^{-2}}{1 + a_1z^{-1} + a_2z^{-2}}$
To this
$H_1(z) = \frac{d_0 + d_1z^{-1} + d_2z^{-2}+d_3z^{-3}+d_4z^{-4}}{1  + c_2z^{-2}+c_4z^{-4}} = H(z)$
since the odd feedback coefficients are 0, it is possible to use pipelining in the feedback loop thus increasing the maximum frequency.
To go from $H(z)$ to $H_1(z)$ ou need to factor
$$H(z) = \frac{b_0 + b_1z^{-1} + b_2z^{-2}}{1 + a_1z^{-1} + a_2z^{-2}} =
\frac{b_0 + b_1z^{-1} + b_2z^{-2}}{1 -2rcos(\theta)z^{-1} + r^2z^{-2}}  $$
$$H(z) = \frac{b_0 + b_1z^{-1} + b_2z^{-2}}{1 -2rcos(\theta)z^{-1} + r^2z^{-2}} * \frac{re^{±(j(\theta+\pi/2))}}{re^{±(j(\theta+\pi/2))}} = \frac{d_0 + d_1z^{-1} + d_2z^{-2}+d_3z^{-3}+d_4z^{-4}}{1  + c_2z^{-2}+c_4z^{-4}}$$
Edit :
If your poles are not complex, the technique I showed will not work.  However you can still apply the look-ahead. You simply need to transform the order-2 filter in 2 cascaded order1- filters and then apply the lookahead transform.
You should consult this site for more information http://people.ece.umn.edu/users/parhi/SLIDES/chap10.pdf


noise remains after applying median filter So I'm apply the median filter to an image, but at the output, there's blue semi-dots appearing. What are they?

This is the output

And the python code for anyone interested

from scipy.ndimage.filters import median_filter
from matplotlib.pylab import imread, imshow, subplot, title, show, get_cmap

img = imread('images/speckle.gif')
new = median_filter(img, 3)

subplot(121); imshow(img, cmap=get_cmap('gray'));
title('The original')
subplot(122); imshow(new, cmap=get_cmap('gray'));
title('Filtered with the median_filter')

show()

Convert the image to gray scale before applying the median filter.


Smoothing a power spectrum My question is regarding how one could/should smooth a rather noisy power spectrum obtained from measurements. There are various methods one can use (a moving average, binning, etc) and I am not sure what the correct approach would be, given the situation I will describe.

Although at first sight it might sound like I am asking about what kind of periodogram to use (see for example this question with excellent answers), I think my situation might be a bit different. What I have is a bunch of complex valued autocorrelation functions (ACF): $n$ datasets of length $m = 8192$, centered at $t = 0$, taken with a sampling rate of 1 nanosecond (so 1 GHz sampling rate).

While the specifics are not too relevant to the question, they can also not be changed; the data is obtained using an FPGA that receives a complex value every nanosecond, and it internally calculates the autocorrelation function using a bunch of steps involving DFT's implemented in the binary hardware.

So the data I have I am essentially stuck with. Given this set of ACFs, I'd now like to get their power spectra, which is of course just the Fourier transform. Here is an example:

There are a few things to see. First, there are two peaks between -100 and 0 MHz. These are the main objects of interest; they change their width and height in each dataset, and I'd like to integrate them to see how their area changes. As they are now they are rather noisy, and I think they would benefit from smoothing.

The second thing is the major peak between 100 and 200 MHz; this peak is always there (always at the exact same frequency with this resolution), and it should be kept. That means that it shouldn't be smoothed, nor should its neighbours take it into account when smoothing. I suppose one solution is just to take it out, smooth, and put it back in.

Finally there is some noise background that varies from about 0 to 0.002. This I am not interested in; it doesn't have to go away completely, but it should still be flattened out so that I can see it as a sort of static background.

So, taking this all into account, I am looking for some advice on what to do. A simple approach would be binning; splitting the data up into $k$ segments of length $l$ and taking the central moment of each segment as my new dataset. Of course it reduces the number of points. Another alternative is using a moving average, which replaces each datapoint by its r-range average. And then there are a million more smoothing methods; I just don't know what is more suitable.

And I suppose this is the main question. Given that I have a dataset that looks like shown above, with the features I described, how should I go about narrowing down what type of smoothing method is most suited? The main deal, as I mentioned before, is finding the area of the two peaks between -100 and 0. Also making the data look a bit more pleasant is not unimportant. The moving average does all this, but it feels so.. primitive.

For reference, the dataset (as comma seperated values) can be found here. Note that this is already the PSD; the Fourier transform of my original datasets.

Two things: 1) Unless there's some kind of synchronicity between your spectrum analyzer oscillator and the signal source, the spectrum should be fairly symmetrical around 0Hz (except for the AWGN). So I strongly suspect something wrong there, because it isn't symmetrical. 2) Since (I assume) you are getting more than one spectral sweep, for each frequency bin take something like the average of the bin over multiple sweeps. For example, a given frequency bin could be the average of the last several sweeps for that bin, or its an IIR of FIR of the bin history. Do NOT smooth across frequency unless you only have one sweep of data and nothing else.


Detecting Pattern from Signal Data by Gaussian Mixture Model? I'm a machine learning newbie.

I have sensor data which is generated by several sensors.

The data is a series of 'time's. (it is not labeled, in other words, I cannot know which sensor generates which 'time').

And each sensor basically generates data periodically, but sometimes it generates other data.

For example, there are two sensors (sensor A generates a 'time' every 2 seconds, sensor B generates a 'time' every 5 seconds.).

Then the data can be: {2.1, 4.3, 5.2, 6.1, 7.9, 9 (noise from one of the sensors), 10.1, 10.2} .... (there are some noises)

What I want to know is that, when another data is given, it is possible to know that the same sensors generate it? (We know the period)

What method should I have to use?

I found the Gaussian Mixture Model. Am I right to use it?

Thanks.

If the data is cyclic by its nature the best thing would work using its spectrum.

You can easily build a system which checks sub set of data to verify periodic and the once you establish your groups checking the affinity of new data is easy - add it to each series. It should belong to the one creates less spread in the frequency (Namely, it follows the cycles).


MFCC feature extraction I have a block of audio and I have performed an FFT on it.  Now what I want to do is convert this FFT into a set of MFCCs, however while I know I need to do something with weighting the FFT samples based upon a set f triangular filters I haven't really got a clue what this involves.

Can anyone go into some good detail on it?  If you are going to give me a matlab example please explain each step well as matlab allows you to cut quite a few corners and I'm trying to implement MFCC extraction on iPhone.

There's a lot of literature on MFCCs on the web, so it would be a bit easier if you could be more specific as to which part of the processing you don't understand. But I'll give an overview of what needs to be done, hoping this is helpful for you:

compute the squared magnitudes of the FFT bins
weigh the bins using triangular windows; usually the windows are chosen such that the centers of the triangles are equidistant on a mel-frequency scale, and such that each triangle begins and ends at the centers of the two adjacent triangles. The mel-frequency scale is defined by
$$m=2595\log_{10}(1+f/700)$$ where $f$ is the frequency in Hz. Look at the figure to see how it works:

Note that the kernels are normalized such that the sum of the weights per triangle equals 1. Usually around 20 such triangular windows are used.
Take the logarithm of the weighted coefficients.
Compute the DCT.


Algorithms to fix audio feedback loop I'm a student and I have been assigned to make a presenation on algorithms to prevent audio feedback loop (microphone feedback, that issue when you hear the voice echoed multiple times).
I looked online but I can't find any specific algoritm name that is used to achive this. I just need a name I can look online for books, whitepapers for that algorithm.
Thank you!

Here is an excellent reference on the topic. It's a survey on acoustic feedback control.
http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.231.9808&rep=rep1&type=pdf


Sampling theory and frequency spacing If a CT signal $x(t)$ is sampled for one ($1$) second and $4096$ samples are generated:

Isn't the highest frequency that can be sampled without aliasing ($4096/2$) $2048\textrm{ Hz}$?
If a $4096$ FFT of this signal is computed what is the frequency spacing between the output points (in $\textrm{Hz}$) ?

Yes, the highest frequency is half the sampling rate. The frequency spacing is 1 Hz since samples are accumulated for 1 second.


Position of poles and Stability in $z$ domain We know in Laplace Transform, if the poles lie on the left of $j\omega$ axis, we can say the system is stable. Similarly can we comment on the stability based on poles position in $\mathcal Z$-transform domain?

Short Answer:
All the poles of a causal (right-sided) and stable LTI system must be inside the unit circle whereas all the poles of an acausal (left-sided) and stable LTI system must be outside the unit circle.
Explanation:
First, consider the following right-sided (causal) exponential sequence $h[n] = a^n u[n]$ whose $\mathcal{Z}$-Transform is:
$$H(z) = \sum_{n=0}^{\infty} a^n z^{-n} = \frac{1}{1 - a z^{-1}} ~~~,~~~ |z| > |a|.$$
As can be seen for a right-sided (causal) sequence the ROC of $H(z)$ extends outward from the pole to the infinity.
Now, consider a causal LTI system with right-sided impulse response $h[n]$ whose transfer function can be expressed in terms of ratio of polynomials:
$$H(z) = \sum_{n=0}^{\infty} h[n]z^{-n}  = \frac{P(z)}{Q(z)} = K \frac{\prod_{k=0}^{M-1} (1 - d_z z^{-1})}{\prod_{k=0}^{N-1} (1 - d_p z^{-1})} $$
Where $d_p$ are the poles of the transfer function at which $H(z)$ diverges. Hence, associated with the $\mathcal{Z}$-Transform $H(z)$, is a region of convergence (ROC) which is singly connected region that does not include any poles inside. Which implies that ROC for the causal system extends outwards from the largest pole.
On the other hand, the BIBO stability of LTI system requires that $|y[n]|<\infty$ for all bounded inputs $|x[n]| < B$, which translates to $h[n]$ be absolutely summable as follows:
$$ |y[n] = \sum_{k} h[k] x[n-k] | < \infty \implies \sum_{k=0}^{\infty} |h[k]| <  \infty$$
This means, for stable system, at $|z| = 1$, the $\mathcal{Z}$ transform
$$ H(z) = \sum_{k=0}^{\infty} h[k] z^{-k} < \sum_{k=0}^{\infty} |h[k] z^{-k}| = \sum_{k=0}^{\infty} |h[k]|  <  \infty$$
converges! Thus for a BIBO stable LTI system, its ROC must include the unit-circle $|z|=1$ inside.
Combinig these two concepts yields: for a causal LTI system to be stable its ROC must include the unit circle and it should extend outwards  from the largest pole to infinity. The corollary of this fact is that the largest pole of the stable and causal LTI system is inside the unit circle. Hence we conclude that all of the poles of the stable and causal systems are inside the unit circle.
An opposite complementary conclusion occurs for the stable but acausal LTI systems (whose $h[n]$ includes left-sided terms like $-a^nu[-n-1]$) for which the ROC is inwards from the smallest pole, and must include the unit circle (for stability), and therefore that all of its poles must be out of the unit circle; i.e., it's smallest pole is larger than $|z|=1$.


Convolution Reverb Calculation Here's my basic understanding of how to implement a convolution reverb:

I measure the IR of a space using a sine sweep. The raw mic output is stored as a Voltage-Time data set. Let's assume that this recording is 5s long.

I run a fft on the recorded mic output to derive the Mag and Phase at every frequency bin. Call this the "IR Buffer".

Next, I have my music source. Let's say it's a vocal track that's playing in real-time. I record this vocal track into a 5s buffer.

Similar to step 2, I run a fft to extract the vocal track buffer's Mag and Phase values at every frequency.

Then, I multiply the Mag and Phase values of the vocal track with that of the IR Buffer at every frequency bin.

Finally, I run a ifft on the multiplied buffer and play the output as the reverb'd vocal track.

Is my understanding above correct?
What I'm very uncertain about is that it looks more like I've created a simple EQ rather reverb.
More specifically, I can't figure out how to extract and apply the time-variant Mag and Phase characteristic of the reverb space. For example, there's a certain decay time for 100Hz in the space, which means that Mag and Phase varies with time for this frequency. How is this time-variant component added to the vocal track?

First, you have to extract the Room Impulse Response (RIR) of the space from your sine sweep recordings. You can do this by taking the ratio of the FFTs of the  measured response and sine sweep input and then taking its IFFT, rir = ifft(fft(measured, Nfft)./fft(sin_sweep,Nfft)). I generally use Nfft = 2^nextpow2(length(measured))
You can truncate the RIR to $L$ samples, depending on how long the reverb tail is.
Break your music source into frames of $M$ samples.
For each frame, you want to convolve a signal of length $M$ with the RIR of length $L$. The output of the resultant convolution will be $M+L-1$ samples. So you have to pick an FFT size accordingly. Then multiply the FFT of the current frame of the signal, with the FFT of the RIR.
The final step is overlap add.

Here is some MATLAB style code.
% sig - Input signal
% rir - Room Impulse Response extracted from sine sweep
% M - length of each frame
% L - length of the RIR

nframes = round(length(sig)/M));
Nfft = 2^nextpow2(L+M-1);
conv_sig = zeros(nframes*M,1);
for frame_no = 1:nframes
    sig_frame = sig((frame_no-1)*M+1:frame_no*M);
    conv_frame = ifft(fft(sig_frame,Nfft).*fft(rir,Nfft));
    conv_frame = conv_frame(1:M+L-1);
    conv_sig((frame_no-1)*M+1:frame_no*M+L-1) = conv_sig((frame_no-1)*M+1:frame_no*M+L-1)+conv_frame
end


Viterbi decoder K=3 using ADD-COMPARE-SELECT Operation I am in trouble with Viterbi decoder for a while. I have some questions.
I need to decode messages from a continuous stream which is encoded.
I'm reading from book : digital communications bernard sklar second edition and the Viterbi decoder including what's called stage "ADD-COMPARE-SELECT " and I didn't understand well how it works , still confused at this stage.
there's branch metric for next possible states and for output table/metric that're related to
ADD-COMPARE-SELECT stage of Viterbi decoder , could anyone please help me out and explain to me in simple words how actually ADD-COMPARE-SELECT stage of Viterbi decoder works? it would be appreciated to explain by actual example.
attachment from the book that Im reading from:
enter link description here
thanks alot.

For each state in the Viterbi decoder, you want to keep track of the shortest path through the trellis that ends in that state. The trellis structure defines what all of the possible state transitions are.
Therefore, each time you get a new symbol observation, you loop through each state in your decoder and you:

Add the metric for each transition that ends in the state to the metric of the previous state, yielding a total metric for each state->state transition,
Compare the list of total metrics, and
Select the minimum metric as the most likely path to the current state.

The Viterbi algorithm, in this way, is really just repeated execution of this add-compare-select process, thus keeping track of the minimum-cost path to each state at any given time.


Would subtracting the phases of two images be a superior difference metric than subtracting the images directly? I'm hoping someone can sanity check this idea as I am admittedly a bit of a noob when it comes to working with FFT's.

Say I have two (registered) images of two of the same object and want to use one as a baseline to check for differences in quality control.  My initial naive approach is to just subtract the two images directly and then whatever remains can be treated as defects.  However, this method is subject to error in the presence of illumination differences.

I'm thinking that I would be better served if I take the FFT of the images, and subtract only the phase information, and then use the IFFT of that result as the defect map, and that this should hopefully eliminate false positives due to lighting.

Does this seem like a reasonable assumption or is there some kind of detail that I'm overlooking?  Thanks for any advice you can offer!

You don't have to work with FFT's and their complex phases, it is far simpler to work directly with the data.  It is particularly simple if the two images are registered.  At the simplest level you find the edges in both of the images and take the difference between the edge maps.  All structure that is common to both images is then zero and all that is left is the pixels that have changed.  This is robust against changes in illumination provided you match the histograms of the edge images prior to subtraction.

Inevitably there are complications, all of which in any case are present when working in the Fourier domain.  I have already touched on the effect of illumination, you can make that slightly more effective by renormalising the histograms of the images before applying the edge filter.  As for which edge filter to use, Roberts, Sobel, or whatever, I would suggest taking the 3x3 Laplacian (this is equivalent to finding the trace of the Hessian of the image density field).

One serious complication is noise: it is essential to remove the noise in the images before starting out.  A simple 3x3 2D smoothing filter passed over the image is fine.  (In the Fourier domain you kill the highest frequencies).  This is important because differentiating random data, which is what edge filters are about, is generally a bad thing to do, and taking second derivatives (the Laplacian) is even worse.  Any noise filter applied in the spatial domain will require a threshold filtering to eliminate noise.  Choosing the threshold visually is effective, but not practical if there is a lot of data.  Then there are auto threshold mechanisms available, and some of those are adaptive to the local noise.  It can get as complicated as you like!

A further complication arises from knowing a priori what is the scale on which you expect to find the difference: how many pixels?  Simple 3x3 filters are inadequate if the features delineated by edges on a scale of >10 pixels.  If you are searching for a particular scale it's worth smoothly shrinking the images so as to make the edges of interest conform to the sought after features.

If you have no idea what the scale is, then my advice would be to use a simple wavelet transform (like the 5-3 wavelet) which automatically locates edges on all scales. I don't have a ready reference to that, but I'll find one if that is of interest.

It can get quite involved depending on how difficult it is to define what you are looking for.  The Multi-scale Morphology Filter was introduced in my field in the paper http://arxiv.org/abs/0705.2072 and refined in a later astronomy paper, http://arxiv.org/abs/1209.2043 .  These dealt with the problem of finding and classifying structure on a point set distribution (galaxies) and, in terms of the question, describe the scale free analysis used to find edges (filaments of galaxies) of arbitrary scale.  The work was also focussed on making the search for structures parameter free, i.e. "objective" and since then has progressed to being able to handle data in N-dimensions of billions of data points (not finished and so not published as yet).

I also use this in industrial process control using thermal imaging cameras, which have weird noise characteristics.

The bottom line is to assess the complexity of your data by doing the very simplest thing contained in the first paragraph.  If that works then don't go any further.


DC Components and Frequency I have started learning digital communication and was on line coding. I read there are certain characteristics that we need to keep in mind while choosing the right line coding method. One of these is considering DC component. I'm not quite able to understand what DC component really is, what and why does it poses a problem in digital communication?

To explain my doubt, I would like to tell two different definitions of DC component (at least they seem different to me) that I read from two different sources. First,

  Signals creates very low frequencies when the voltage level in a digital signal is constant for a while. These frequencies around zero are called DC components, creates problems for a system that cannot pass low frequencies or system that uses electrical coupling.

Second,

  Over one cycle period of a waveform, if all the positive voltages are cancelled by negative voltages, then, DC component of the waveform is 0. In line coding, the signal with the non-zero DC component is treated as a distorted one and it can create errors in received signals.

I don't seem to have any clue about what the first definition wants to say. As far as I know, the definition of frequency is number of cycles completed per unit time. Actually, how do you even measure the frequency of a digital signal? And how does it relate to the DC component?

I understand the second definition in itself pretty well but I don't understand how does it relates to the first definition. Also, what problem does it poses in transmission if the positive voltages do not cancel the negative ones?

P.S : A less mathematical reply would be highly appreciated since I'm new to this field.

Thanks in advance.

how do you even measure the frequency of a digital signal?

This may be the source of your confusion. In digital communications, data is transmitted using an analog, continuous-time waveform. This is necessary because only this kind of waveform can exist in the physical world. Like all waveforms, it has a bandwidth, and its Fourier Transform determines what frequencies are present in it.

In digital communications, the waveform used to communicate has a specific form: it is a sequence of pulses, and the amplitude of the pulses encodes the transmitted bits.

The simplest waveform is one that encodes a bit $1$ with a pulse of positive amplitude, and a bit $0$ with a pulse of equal but negative amplitude. You may think that, for random bits, the negative pulses cancel out with the positive pulses and the average DC level of the signal is 0. This is in general not true over short time periods; even random bits will have sequences with many more $1$'s than $0$'s, or viceversa. This differential causes the signal's DC to grow, which causes distortion at the receiver end.

Line encoding can help with this problem. For example, the Alternate Mark Inversion line code encodes a bit $0$ to a pulse of amplitude zero, and a bit $1$ to pulses of positive and negative amplitude, alternating. So, every positive pulse is followed by a negative pulse, and the signal's DC never builds up.


Calculating Power of Sinusoidal Term and a Convolution expression The following is an excerpt from the book Statistical Digital Signal Processing and Modelling (Monson Hayes)

  ...consider a random process consisting of a random phase sinusoid in
  white noise $$ x(n) = A\sin(n\omega_{0} + \phi) + v(n)$$ where $\phi $
  is a random variable that is uniformly distributed over the interval
  $[-\pi, \pi]$, and $v(n)$ is white noise with a variance of
  $\sigma_{v}^{2}$. The power spectrum of $x(n)$ is

  $$ P_{x}(e^{j\omega}) =  \sigma_{v}^{2} + \frac{1}{2}\pi
 A^{2}[u_{0}(\omega - \omega_{0}) + u_{0}(\omega + \omega_{0})]$$

  Therefore it follows from Eq. (8.23) that the expected value of the
  periodogram is  $$E\{{P_{per}(e^{j\omega})}\} =
 \frac{1}{2\pi}P_{x}(e^{j\omega})*W_{B}(e^{j\omega})$$ $$=
 \sigma_{v}^{2} + \frac{1}{4}A^{2}[W_{B}(e^{j(\omega - \omega_{0})}) +
 W_{B}(e^{j(\omega + \omega_{0})})]$$

I am unable to figure out how he obtained the expression for power of the sinusoidal term. Shouldn't it be just $\pi^{2}A^{2}$ multiplied by the bracket term?
Considering that $W_{B}$ represents fourier transform of Bartlett window, how is the convolution expression being evaluated here?

The power spectrum $P(e^{j\omega})$ is the Fourier transform of the autocorrelation function of $x(n)$:

$$R_x(k)=E\{x(n)x(n+k)\}\tag{1}$$

If you evaluate $(1)$ you should get

$$R_x(k)=\sigma_v^2\delta(k)+\frac{A^2}{2}\cos(k\omega_0)\tag{2}$$

Taking the Fourier transform of $(2)$ leads to the given expression for $P(e^{j\omega})$.

The convolution is obtained by noting that

$$u_0(\omega-\omega_0)\star W_B(e^{j\omega})=W_B(e^{j(\omega-\omega_0)})\tag{3}$$

and

$$\frac{1}{2\pi}\sigma_v^2\star W_B(e^{j\omega})=\sigma_v^2\frac{1}{2\pi}\int_{-\pi}^{\pi}W_B(e^{j\omega})d\omega=\sigma_v^2\tag{4}$$


Stereo triangulation similar triangle derivation I understand that the similar triangles are the red and blue triangles. I'm trying to understand where $T + x_r - x_l$ comes from. I understand we are trying to compute the length of the base of the blue triangle, but I'm not sure how this expression is finding that.

You are computing the base of the red triangle, not the blue one. The blue triangle's base is already indicated as $T$.

Not very clear from the given figure, but as the only consistent answer I could find, assume that the coordinate of the right corner, $P_r$ of the blue triangle $P_l P P_r$ is measured with respect to the second origin $O_r$ indicated in the bottom right, and it seems to be a negative number by definition.

Therefore the length of the base of the triangle $P_l P P_r$ can be given as
$$ L = |T| - |x_l| - |x_r|$$

where $T$ is a positive quantity by default, $x_l$ is measured wrtt left point $O_l$ and is positive by default, and $x_r$ will be a negative number as stated above; hence the absolute values can be expanded as

$$ L = T - x_l + x_r = T + x_r - x_l$$

Note that actually you do not have to use absolute values, but just proper adition of the signed quantities would suffice to derive the result. I've used the absolute values to provide some intuition into the signs of the summing coordinates.


Differences Using Maximum Likelihood or Maximum a Posteriori for Deconvolution / Deblur Are there any differences if you use Maximum Likelihood or Maximum a Posteriori to estimate the Point Spread Function for image deconvolution?

You can think of MAP as a regularization of the ML.
Just like you have regularization for Least Squares Problem (They can be built, mostly, as MAP problem).

The nice thing is that, as always, the best regularization is more data, namely, in most case when there is a lot of data they collide (Namely, low sensitivity fir the Posterior PDF).

So they differ mainly when there is (Relatively) low amount of data.
Now, when you have low amount of data and you know nothing about the parameters you're trying to estimate, ML is the way to go.
If you have low amount of data yet some prior knowledge about it, or reasonable assumption, make those assumption as regularization.
It even good to make the effort and describe this knowledge in Posterior PDF form.

Update

Found really nice tutorial about the topic - The Truth About Priors and Overfitting.


How to 'interpret' the Fourier Transform (specifically, of a convolution kernel) As part of a homework assignment, I had to take the Fourier transform of the kernel I was using to convolve a signal.  The kernel was a constant rectangular function, that was 1 within the square $(-1, -1)$, $(1, 1)$ and 0 everywhere else.

I was wondering what this result really means.  What have I gotten once I take the Fourier transform of this function?

Convolving a signal with a kernel is equivalent to applying a filter to it. The kernel is the impulse response of the filter, and the Fourier transform of the kernel is thus the frequency response of the filter.

In your case, the filter's impulse response is a rectangular function of width 2 and centered at 0. You can interpret this filter this way: for any $t$, the filter output $y(t)$ will be the average of the input signal $x(t)$ over the interval $[t - 1, t + 1]$. This is an averaging, low-pass filter, and indeed, the Fourier transform of the kernel is a decreasing function, showing that the higher frequencies are attenuated.


Digital Signal Processing Filters Problem Hello I need some help with this question.
$$[] = _1[ − _1] + _2[− _2] + _1 [− _3]$$
Set $a_1, a_2$, and $b_1$ to any positive or negative real numbers of your choice. Set $k_1, k_2$, and $k_3$ to any positive integer numbers of your choice ($k_1$ can be zero). Present this equation with the numbers you have chosen. Present the transfer function $()$ of your filter
I have chosen the values so my equation looks like this
$$y[] = 0.4[] + 0.6[− 1] + 0.8[− 2]$$
I have worked out $H(z)$ to be $$\frac{0.4 + 0.6z^{-1}}{1-0.8z^{-2}}$$   is this correct?
My questions are:

Is this filter an FIR or IIR type?
How do you determine the poles and zero?
How do you determine if the filter is stable or unstable?

Thanks!

I have worked out H(z) to be (0.4 + 0.6z^-1) / (1-0.8z^-2) is this correct?

Close but one detail is wrong

Is this filter an FIR or IIR type?

Look at the pole locations

How do you determine the poles and zero?

Poles are the roots of denominator polynomial. Zeros are the roots of the numerator polynomial.

How do you determine if the filter is stable or unstable?

Look at the pole locations
These questions are about as fundamental as it gets in DSP. They are also incredibly simple, in fact, so simple that it's hard to give a hint. If it takes you more than a few seconds to answer this, you need to adjust the way your are studying. I strongly recommend to not proceed to any more advanced topic until you very comfortable solving this type of problems. It's fundamental to whatever comes next.


How to Add Group Delay? I was confused about the value of group delay, I am not sure whether my opinion is correct. For example, if I want to add a group delay of $0.25$ time units to a zero-phase signal, does it mean that I just need to multiply the signal by $0.25\omega$ in frequency domain?

This is a part of my lab of DSP. In this lab, I was asked to design a LPF by least squares method. In the first question, $h_d(\omega)=1$, and it is asked to have a $0.25$ time units group delay in the next question. I tried it but I found nothing different happened, so confused.

You got the shift a bit wrong – a shift in time domain (which you want) is a multiplication with a sinusoid in frequency domain.

This is simply a consequence of the convolution theorem: you want to convolve with a dirac $\delta(t-\Delta t)$ in time, so you need to multiply with the dirac's Fourier Transform in frequency domain.

Thus, what you'd need is

DFT -> $\cdot e^{j2\pi \tilde{\Delta t} f}$ -> IDFT.

Yet, you incur other problems here:

To do a DFT of your whole data, you'd need to first put your whole signal in a buffer. In reality, that means you can't do that to work with live signals.
The DFT will lead to a cyclic temporal shift. I.e. the (fractional) samples that you "shift out" at the end will appear at the start of your output! That means you can't simply divide your signal into smaller DFT chunks and shift each one of them separately.

If you need to avoid that: Use the fast convolution; overlap-add and overlap-save are the two algorithms you should be looking at. It's what I did in a lot of places, when I needed multiple variable delays.

Generally, any (time-)symmetric linear filter has constant group delay – you don't even have to go through frequency domain if you don't like to.


About frequency resolution The frequency resolution of DFT is

$\Delta f = \frac{1}{T_{0}} = \frac{1}{NT} = \frac{f_{s}}{N}$

where $f_{s}$ is sampling frequency, $T_{0}$ is sampling time, $N$ is number samples, and $T$ is the sampling interval.

Where does this definition come from? Espescially the $1/T_{0}$, why this equality?

The longer the total sample window (T0), the smaller the frequency difference between two sinusoids with exactly M periods and exactly M+1 periods (for integer M <= N/2) in that window, which is required for orthogonality between adjacent DFT bins.

Also, the bin spacing is not the resolution. The frequency resolution can be greater or less than the DFT bin spacing, depending on the S/N ratio and your resolution needs (separation or estimation).


Equal power crossfade I am working on an audio system where I need to implement an equal power crossfade between two signals over a specified number of samples. Initially I implemented a linear fade but this creates a dip in the middle of the crossfade which I don't want. I tried a few other methods that I found on Google but they either make a dip in another part of the curve or even make the sound louder.

If I have float samples ranging from -1.0 to 1.0 and I want to crossfade between them over the course of x samples, keeping the perceived volume the same over the entire crossfade, what would be the correct formula for calculating the curve? Preferably I want to have a value like 0.001 represent the minimum level that the fade in starts at and the fade out goes to. 1.0 should be the maximum. Thanks in advance for any advice!

if the two audio signals are totally uncorrelated, the the squares of the two crossfade gains should add to 1.  so, to prevent that dip, your crossfade function will look more like $\sqrt{\frac{1}{2}(1+t)}$ and $\sqrt{\frac{1}{2}(1-t)}$ rather than $\frac{1}{2}(1+t)$ and $\frac{1}{2}(1-t)$  for $-1 \le t \le +1$.  note, for constant power crossfades, at the midpoint of the crossfade ($t=0$), both crossfade gains are at $\sqrt{\frac{1}{2}} \approx 0.707$ rather than at $\frac{1}{2}$.  that's what addresses the dip you want to avoid.

if your two audio signals are somewhere between totally uncorrelated and totally correlated, you need to determine the degree of correlation (that's what cross-correlation is all about), and then the crossfade functions can be determinied analytically.  i have posted to music-dsp a little theory about it.  it's not perfectly rigorous, but i think it works.


Inverse DFT in terms of original signals for product of two DFTs Given $x[n]$ and $h[n]$ that are zero for $n<0$ and $n>15$, their linear convolution being $y[n]=x[n]*h[n]$ and their discrete time Fourier transforms are $X(e^{j\omega})$ and $H(e^{j\omega})$, then
$Z[k]=X(e^{j\omega})H(e^{j\omega})|_{\omega=\frac{4\pi k}{31}}, k=[0,...,30]$.

How can $z[n]$ be expressed in terms of $y[n]$, where $z[n]$ is 31-point inverse DFT of $Z[k]$.

Again in this problem you can make use of the duality property to derive the output result without much involvement in the DFT sum formula. Before beginning, however, it's necesary to remember the time-frequency relations of two basic operations, namely, the expansion and compression blocks.

The EXPANSION block is the following:

$$ x_N[n] \longrightarrow \boxed{ \uparrow 2 } \longrightarrow y_{2N}[n] $$

shown operating on a time-domain signal $x_N[n]$ of length $N$ and producing the output $y_{2N}[n]$of length $2N$ given by the relation :

$$y_{2N}[n] = \begin{cases}{ x_N[n/2] ~~,~~~n=2m ~,~~~~~~ m=0,1,..,N-1 \\ ~~~~~0~~~~~~~~,~~~ n = 2m+1 ~, m=0,1,..,N-1} \end{cases} \tag{1}$$

The frequency domain relation between $N$-point DFT $X_n[k]$ of $x_N[n]$ and $2N$ point DFT $Y_{2N}[k]$ of $y_{2N}[n]$ is:

$$ Y_{2N}[k] = X_N[ (k)_N] ~~~,~~~ k = 0,1,2,...,2N-1 \tag{2}$$

Verbally stated, $Y_{2N}[k]$ is a cascade (mere repetition) of $X_N[k]$ and can be shown with the following matrix symbolism, each X,Y are (say) row vectors :

$$ Y = [ ~~~-~ X ~ -~~ | ~~~-~  X ~~~-~ ] \tag{3}$$

The COMPRESSION block is the following:

$$ x_N[n] \longrightarrow \boxed{ \downarrow 2 } \longrightarrow y_{N/2}[n] $$

shown operating on a time-domain signal $x_N[n]$ of length $N$ and producing the output $y_{N/2}[n]$of length $N/2$ given by the relation :

$$y_{2N}[n] = x_N[2n] ~~~, ~~~ n = 0,1,..,N/2-1 \tag{4}$$

The frequency domain relation between $N$-point DFT $X_N[k]$ of $x_N[n]$ and $N/2$ point DFT $Y_{N/2}[k]$ of $y_{N/2}[n]$ is:

$$ Y_{N/2}[k] = 0.5( X_N[ k ] + X_N[k + N/2]) ~~~,~~~ k = 0,1,2,...,N/2-1 \tag{5}$$

Verbally stated, $Y_{N/2}[k]$ is obtained from $X_N[k]$ by folding: adding second half of $X_N[k]$ onto its first half.

$$ Y = [ ~-~ {}_{N/2}X ~-~] + [~-~ X_{N/2} ~-~ ] \tag{6}$$

Then the DUALITY principle states that, if frequency domain sequences are processed by those expansion and compression blocks, then their effects in the time-domain seuqences will be the same as above that happened on the frequency-domain sequences of time-domain expanion and compression. Note that, when different sample rates (or number of DFT points) are used in the input and output of the processing blocks, which is the case for the expander and compressor examples, then the duality relation is scaled by the ratio of sampling rates which would cancel the 0.5 at the Eq(5) in our case.

By using these concepts, we can describe the time-frquency domain relation of the problem you have provided.

Now observe that your output signal $$Z[k]=Y(e^{j\omega})|_{\omega = \frac{2\pi}{N} 2k} = [Y(0) , Y(2), ... , Y(30) | Y(1) , Y(3) , ..., Y(29) ]$$ can be given by the following block diagram (with $N=31$):

$$ Y_N[k] \rightarrow \boxed{ Y_N -|- Y_N } \rightarrow W_{2N}[k] \rightarrow \boxed{ \downarrow 2} \rightarrow Z_N[k] $$

Notice that the first blcok is given in (3) and related to the expander operation, where as the second block is clearly a compressor. Lets denote the time domain signals by $y_N[n]$, $w_{2N}[n]$ and $z_N[n]$, and proceed from left to right as follows:

First, given the relation beween the DFT sequences $Y_N[k]$ and $W_{2N}[k]$ then their time-domain originators is the expander operation and hence we have from Eq(1):

$$w_{2N}[n] = \begin{cases}{ y_N[n/2] ~~,~~~n=2m ~,~~~~~~ m=0,1,..,N-1 \\ ~~~~~0~~~~~~~~,~~~ n = 2m+1 ~, m=0,1,..,N-1} \end{cases} \tag{7}$$

And secondly, given the compressor block applied on the frequency signal $W_{2N}[k]$, we conclude its effect, by using duality, on the time domain signal $w_{2N}[n]$ as:

$$ z_{N}[n] = w_{2N}[n] + w_{2N}[n+N/2] ~~~,~~~n=0,1,...,N-1 \tag{8}$$ namely the folding operation on $w_2[n]$.

Combining Eq(8) and Eq(9) we conclude the relation between N-point sequences $z[n]$ and $y[n]$ as:

$$\boxed{z[n] = \begin{cases}{ ~~~y_N[n/2] ~~~,~~~n=2m ~,~~~ m=0,1,..,N-1 \\ y_N[n/2 + N/2],~~~ n = 2m+1 } \end{cases} }$$

Eventhough the description seems really long and somewhat abstract, the operation is really neat and quick as the following Matlab/OCtave code demonstrates  it:

N = 31;
y = 1:N;   % just some signal

Y = fft(y,N);    % obtain YN[k]
W = [Y, Y];      % W = [ Y | Y] , w[n] = y[n/2]: is expanded verison of y[n]
Z = W(1:2:end);  % Z[k] = W[2k]; compress W[k] to get Z[k]

z = ifft( Z, N); % obtain z[n] from Z[k]

% Part-II: Use the time-domain relation as described in the answer to obtain z[n]
we(1:2:2*N) = y ;           % we[n] = y[n/2]: expanded version of y[n]
we = [we , 0];              % Append (2-1) zeros to the end of we if N is odd.
z2 = we(1:N) + we(N+1:2*N); % folding we[n] to obtain we[n]

% Display:
[z' , z2' ]


Can I measure non-linearity of a line echo path? I'm dealing with an echo signal from a hybrid. An echo path includes DAC, the hybrid, ADC and the line infrastructure. I can see, that after using a digital linear echo-canceller the residual echo signal is still significant and I suspect that echo-path has additional non-linear components. Now I want to estimate, how much of the incoming signal is an actual echo and not a noise, and I want to do it before actually implementing or inferring any specific non-linear function (e.g. without assumptions regarding the order of Volterra series).
I've found the ITU-T O.42 recommendation, which suggests using the 4-tone intermodulation method to measure nonlinear distortion impairments on analogue circuits. But it doesn't suit me, since my echo path includes the DAC.
So is there a way to measure or estimate the level of the non-linearity in my echo path? I would highly appreciate any references. Also, if this question is more related to electronics, please feel free to migrate it to the corresponding site.
UPDATE:
Definition of non-linearity: based on the source, a non-linearity for an echo signal is anything that cannot be modeled with a linear filter: a dc component, a long tail and truly non-linear dependencies. So the question becomes: is there a way to measure or estimate the level of those components compared to noise? The peculiarity here is that both a dc component and a long tail are actually linear in their nature.
Signal measurements: Following are the data that I got by changing the power of transmitted signal (from 1 to 6).
Power spectral density plots (regular and normalized), where frequencies are limited by a sampling rate:

Left: echo mean z mean and variance z var and residual echo mean e mean and variance e var for a regular (random) transmitted signal (using some arbitrary length of a linear filter). Right: echo mean and variance and residual echo mean and variance (using filter weights from training on the regular signal) for a constant  transmitted symbol.

Observations:

The PSD plots look identical, except for unclear shift of a plot for the lowest power level (1).

Mean and variance plots do not resemble any obvious non-linearity, which would have shown as a power law dependency everywhere. Based on increased e mean for a constant transmitted symbol, there is a tail in the echo.

What puzzles me the most is a linear growth of a residual echo variance, which I cannot explain.

So the original question is still relevant. Additionally, can I actually use obtained plots to estimate the echo non-linearity and if so, how should I do it?

excellent, you have PAM-4!
Send a long sequence of your least-amplitude symbol, and integrate the received power to calculate the energy of your received signal.
Do the same with the three other symbols.
Now, you can solve the equation
$$E_{RX} = \alpha E_{TX} + E_N + F(P_{TX}),$$
with $\alpha$ being the channel gain, which has to be the same for all 4 of your measurements, $E_{RX},E_{TX}$ being the transmitted and received energies (hint: normalize them to the least-power symbol), $E_N$ is noise energy, i.e. noise power times number of samples (which becomes asymptotically constant if you just average enough symbols), and $F$ is the nonlinear error term.
With $E$ known, and $\alpha$ and $E_N$ overdetermined determined from four measurements, $F$ is the reason for the error.


Segmentation of overlapping noisy parallelograms in a binary image Edited.

Previously this question was named "Segmentation of overlapping pulses in a binary spectrogram" and it was stated in slightly different way. I'm apologise to those users, who were affected by reformulation of the question.

Problem.

There are a set of binary images, that contain arbitrary set of noisy parallelograms. The height and width of objects can vary in a wide range. Some parallelograms can overlap. Set of typical images is listed bellow.

I need an algorithm to segment objects on those images in order to prepare them to extraction of connected components. There is one restriction on the segmentation algorithm - shape of the objects should be preserved during segmentation.

Attempted solutions.

Morphological erosion.
Morphological erosion don't work well in this case. In order to separate three overlapping rectangles in the last image from typical set, I needed structuring element so big, that it erases some small objects from other images.

Boundary extraction and Hit-or-miss transform with corner-like structuring elements.
I have tried detect corners of parallelograms, in order to find regions of intersection. But pattern matching of boundary with corner-like structuring elements turned out to be a bad idea (boundary doesn't have the shape of a corner in regions of intersection).

Some other morphological transformations.
I have tried playing with combinations of original and shifted images.
$A - [A - [(A)_{(-x, 0)}\cap(A)_{(+x, 0)}]]$, where $A$ is an original image and $x$ is a shift along $x$-axis. This transformation can be applied for separation of objects on the last image.
But it messes up second image.

Question.

Is there good general approach to such problem?
Any ideas, hints, links, articles will be much appreciated.
Thank you.

After couple days of research, I have found another interesting solution. The Similar algorithm for segmentation of overlapping objects is described in this article (the algorithm below is an adaptation for objects of parallelogram shape). It uses properties of the boundary to perform separation. Perhaps someone may find it interesting.

Overlapping rectangles.

To get the intuition behind the algorithm, consider easy case of two overlapping rectangles.

We can extract the boundary ( by applying morphological transformation $A - A\ominus B$) and choose an arbitrary point on it.

Now we can traverse boundary starting with chosen point and plot the dependence $\lvert \frac{dy}{dl}(l)\rvert$ along the boundary. (Here $y$ is $y$-coordinate of a pixel and $l$ is a length of the path along boundary). Each peak on the plot is a vertical edge of the rectangle. We can select begin and end points of edges by finding intersection of the graph with threshold $t = 0.5 \underset{l}{\max}\lvert\frac{dy}{dl}(l)\rvert$.

After all, we can connect each pair of points with the amount of points on the path between them greater than 3, and distance between them $d < d_{0}$.

More complex parallelograms.

This algorithm is also suited for handling more complex overlapping parallelograms. In this case dependence $\lvert\frac{dy}{dl}(l)\rvert$ will have the same form as for case with rectangels, but reversed around $x$-axis and shifted. So, the same technique can be applied to extract corner points.

Noisy rectangles.

For noisy objects (as in my case) the dependence $\lvert\frac{dy}{dl}(l)\rvert$ doesn't have such straight shape, as in previous cases. So, additional effort can be required. In order to increase performance of an algorithm, I have applied closing to the original image and moving average filter to the sequence $y[l]$, before calculating derivative.


Jack buffer size and FFT frequency resolution In my actual Jack Audio setup (almost default values), Fs = 48kHz the buffer size is 1024 samples.

I'm planning to do some filtering (lowpass, highpas, ...) and choose to use FFT.

Had done some reading and found that the frequency resolution is $ Fs \over n $ so in that case the frequency resolution will be

$$
{Fs \over n} = {48000 \over 1024} = 46,875
$$

Which is a high value and well imprecise value do deal with it.

What if want do 60Hz highpass filtering ? (c[0][1] = 46Hz, c[0][2] = 92Hz c is complex returned from fft function.)

Some details:

I have choose FFT because i think more simple to do such simple things, but i'm a really beginner in the DSP field this can be false.

Ps: In jack you register a callback and it gets called when there's data to be processed, with a buffer and it's buffer size (fixed).

As a general rule which is easy to remember: The frequency-domain resolution of a FFT does only depend on the time-duration of your signal (this rule also works when you do zero-padding of your signal). The resolution is always given by $1/T$, where $T$ is the duration of your FFT window in seconds. So, in your case your frequency domain resolution will be $1/(1024/48000)=46.875\text{Hz}$ as you have pointed out. Hence, if you want to get better frequency resolution with an FFT, you need to have longer FFT windows. You can use zero-padding to reduce the frequency bin distance, but this will not buy you more frequency resolution. You can have a look at this article I wrote about FFT resolutions, zero-padding and other effects.

However, note that you should not use the FFT to perform filtering of a streaming signal: reference 1 reference 2. Instead, you should perform time-domain filtering with a FIR or IIR lowpass filter. Then you have less delay and would also not need to concern about frequency domain resolution. Let me know, if you need an example on how to design such filter in Python.


How to find phase difference between 2 signals, one has an unknown function I have a system that is forced by a sinusoidal function with specific amplitude, phase and frequency ( all these are known ahead of time). I get a deformation of sinusoidal form as well (think of a damped harmonic response), and I would like to find the phase lag. I was advised to use correlation but I don't have the function that produces the response, only gathered data. This what I am doing but honestly I am not confident in my findings:

    y =[] %this is the gathered data, first column is time, second is deformation
    frequency =1000; %Hz,i know this for a fact
    t=y(:,1);        %just getting the time to produce a sine wave for calibrating

    fs =  14*frequency; %this is random, i m just making sure i m above 2*frequency.
    Nfft = fs ;

    sine_wave = 1+sin(2*pi*frequency*t + pi/2);%I need to 1 + for physical reasons

    FFT_sine_wave =fft(sine_wave,Nfft)/Nfft; %FFT of sine_wave
    FFT_Y2 =fft(y(:,2),Nfft)/Nfft;           %FFT of response

    magsine_wave = abs(FFT_sine_wave);       %magnitude of sine_wave
    magY2 = abs(FFT_Y2);                     %magnitude of response

    f=linspace(-fs/2, fs/2, Nfft);           %frequency vector

   [old_max_sine_wave index_sine_wave] = max(magsine_wave);
   magsine_wave(index_sine_wave) = -Inf;  %getting rid of the DC offset

    [old_max_Y2 index_Y2] = max(magY2);
    magY2(index_Y2) = -Inf;

    [maxAmp_sine_wave index_sine_wave] = max(magsine_wave); %Finding new amplitudes
    [maxAmp_Y2 index_Y2] = max(magY2);

    new_ang = angle(FFT_sine_wave);    %getting the angles of sin wave

    ang_miew = new_ang(8002)*(180/pi)+180 % 8002 is where f = 1000 Hz falls in frequency
                                          %vector, so i use the index to find angle

    amp_miew = 2*abs(FFT_sine_wave(index_sine_wave)) %this gives me the correct amplitude

     new_angY2 = angle(FFT_Y2)*(180/pi);  %same reasoning as above
     ang_miewY2 = new_angY2(8002)+180

Please help me out, if any of you need the y vector that gives the time and deformation please let me know. I am in dire need of help!
Thank you in advance

I don't fully understand your question, or what your code is supposed to do. However, this line is wrong:

fs =  14*frequency;

You should use the same sampling frequency as was used to acquire the data. Assuming uniform sampling, you can do something like

dt = t(2) - t(1);
fs = 1/dt;

Also, be aware that estimating the actual phase from the FFT is not straigthforward. Here is a quick and dirty example showing one way to do it. Start with this code:

fs = 1000;
t = 0:1/fs:0.1;
f = 100;    % signal frequency
phase = 0;  % signal phase
y = cos(2*pi*f*t+phase);
Y = fftshift(fft(y,2^nextpow2(length(y))));
Ym = abs(Y);
[m, idx] = max(Ym);
Ym(idx) = 0;
[m, idx] = max(Ym);  % index of f
Ya = angle(Y);
ph_e = Ya(idx);
fprintf('Actual phase is: %1.5f\n', phase);
fprintf('Estimated phase is: %1.5f\n', ph_e);

If you run this code you'll get:

Actual phase is: 0.00000
Estimated phase is: -0.49087

The phase estimate depends on many factors, including the start and end times of your time vector. The thing is that, in this particular setup, the FFT interprets a phase of zero as -0.49087. Then, modify the third to last line in the script:

ph_e = Ya(idx) + 0.49087;

Now the script will properly identify the phase. Try it by setting the phase in line 4 to, say, pi/3, then run the script:

Actual phase is: 1.04720
Estimated phase is: 1.03347

Which is pretty good accuracy.

I hope this helps you to get going.


Obtaining signal energy modulation for voice activity detection Im currently working on VAD (voice activity detection) topic. Since I'm beginner in digital signal processing, I've started to wonder if calculating single frame energy (after aplying FFT to it) and looking for values around 4Hz and bandwith (for example) from 200Hz to 600Hz would be as effective as checking for 4Hz modulation energy after finding MFCC coefficients?

From my understanding of your question, you want to ask whether one should look for energy(of FFT coefficients) of bands directly, or one should use MFCC's (energy) for  making VAD decisions.

As MFCC based approach is to mimic the human auditory systems response, it will be better to use MFCC coefficients.

if one uses same bands as MEL Filterbank for calculating band-wise energy (after taking FFT), than the results will be similar to using MFCC'S, as MFCC's are just the transformation from FFT coefficients.


What is the application difference between extent and solidity in image processing? I've come across the image object properties, "extent" and "solidity."

Definitions are pretty straightforward to understand:
Extent of an image object is defined as area of the image object divided by the area of its bounding rectangle.
Solidity (convexity) of an image object is, area of the image object divided by area of its convex hull.

However, coming to think of these two properties in terms of the applicability, I can't seem to think how one would give an edge compared to the other. After a quick look-up of these terms, my thoughts are that solidity - as compared to extent - can be a good criteria of image object's roundness since the convex hull forms the minimum enclosing polygon around the object.

What other applications do these two terms have - in conjunction with each other or as building blocks for more advanced contour properties?

I would suggest that Solidity is the better measure. Extent would be a cheap approximation. It’s cheap because computing the bounding box is cheaper than the convex hull (and really simple to implement).

Solidity is sometimes referred to as Convexity, but more commonly Convexity is the ratio of the perimeter of the convex hull to the perimeter of the object. These two measures are of course related, but not the same.

Solidity doesn’t quantify roundness. A rectangle has a Solidity of 1 (the maximum posible value). And so does a triangle. Roundness (or similarity to a circle) can be conveniently computed by comparing the area and the square of the perimeter, or the area to the product of the width and height (look for Feret diameters to compute these). Another approach is to determine the coefficient of variation of the distance of each boundary pixel to the geometric center of the object.

Solidity is useful to quantify the amount and size of concavities in an object boundary. Holes are also often included. For example, it distinguishes a star from a circle, but doesn’t distinguish a triangle from a circle.

Here is a good survey of object shape features: M. Yang, K. Kpalma and J. Ronsin, "A Survey of Shape Feature Extraction Techniques", in: Pattern Recognition Techniques, Technology and Applications, P.Y. Yin (Editor), I-Tech, 2008.

The DIPlib library implements the more meaningful of these features. (Disclosure: I implemented most of those features in DIPlib.)


Sampling rate of a sinusoidal Signal What is the minimum sampling rate of a $\sin(\pi t)$ signal? When the sampling frequncy is twice that of the frequency, the spectrum is not present and in time domain all the samples we get is 0. In this case does the sampling rate should be greater than $2 f_\text{s}$? Then how does the Nyquist rate fit in this?

What is the minimum sampling rate of a $\sin(\pi t)$ signal?

$$F_\text{s} = 1 + \epsilon$$

where $\epsilon > 0$.  you must sample at some frequency exceeding twice the bandwidth.  (all this assmuming real $x(t)$ and none of this bandpass sampling stuff.)

  When the sampling frequncy is twice that of the frequency, the spectrum is not present and in time domain all the samples we get is 0.

nooooooo! that's not the case!  you could miss-sample it by $90^\circ $ and

$$x[n] = x(nT) = \cos(\pi n) = (-1)^n$$

instead.

  In this case does the sampling rate should be greater than $2 f_\text{s}$? then how does Nyquist rate fit in this

two things:

the original Nyquist paper (i'll try to find a link to it), actually excludes pure sinusoids, particularly a pure sinusoid sitting at Nyquist.
to deal with the $x[n] = x(nT) = \cos(\pi n + \phi)$ thing (a sinusoid at the Nyquist frequency), you have to accept the notion of Dirac delta functions, $\delta(f)$ in the frequency domain.  and then you will see that, at Nyquist, ($f = \frac{F_\text{s}}{2}$), there will be overlapping and aliasing.  but unlike other aliasing, you lose only half of your information, just the phase.  so if there were reconstruction done, it would be only on the basis of the quadrature component of the

$$ x(t) = cos(\pi t + \phi) = \cos(\phi)\cos(\pi t) - \sin(\phi)\sin(\pi t) $$

input. if $F_\text{s} \triangleq \frac{1}{T} = 1$ and $x[n]=x(nT)$, only the $\cos(\phi)\cos(\pi t)$ term counts.


How to do base-band simulation of coherent detection of M-ary orthogonal signals? I currently have a complex base-band simulator for non-coherent detection of M-ary orthogonal waveforms (such as M-FSK).  I run the noisy signal through M correlators at lag=0, compute the magnitude squared of the results and pick the largest.  Everything is perfectly synchronized.

It is just: modulate -> add AWGN -> demod.

Currently, my non-coherent simulation BER performance matches the theoretical equations derived in Proakis.

What do I have to tweak in order to simulate coherent detection given that my received signal is perfectly synchronized?

Since the answer was given by vaz in the comments I figured I would just post it as an actual answer...

To perform coherent detection for M-ary orthogonal signals in a complex base-band simulation simply take the real part of each matched filter (correlator) output and select the symbol with the largest value.

Also, don't make my mistake and use the absolute value of the real part.  I made those changes and verified against the theoretical coherent equations MATLAB uses.


Reducing Pepper Noise from an Image For an image with pepper noise, presuming it is desirable to reduce the pepper noise, which operations, in which sequence, would you apply? Find an image, add pepper noise, and demonstrate your results.

Salt and Paper noise are better dealt and reduced using Median Filter.

The properties of the Salt and Pepper noise make it an outlier in almost any patch of the image.
Hence the best way to remove it is using a method robust to outliers.
Linear filters (Low Pass Filters) are basically weighted mean an hence very sensitive to Outliers which makes them less efficient in this case.

The Median is known and simple robust estimator and hence perform very well in this case.

The way to use it is simple, working by patches / windows.
Iterate through the pixels of the image, for each pixel open a Window of its neighborhood controlled by the parameter which is the Window Radius.
For each window replace the value of the center pixel by the median of the window.

MATLAB Code (From StackExchange GitHub):

% Image Salt and Pepper Noise Filtering
% See http://dsp.stackexchange.com/questions/27298

%% General Parameters and Initialization

clear();
close('all');

% set(0, 'DefaultFigureWindowStyle', 'docked');
defaultLooseInset = get(0, 'DefaultAxesLooseInset');
set(0, 'DefaultAxesLooseInset', [0.05, 0.05, 0.05, 0.05]);

titleFontSize   = 14;
axisFotnSize    = 12;
stringFontSize  = 12;

thinLineWidth   = 2;
normalLineWidth = 3;
thickLineWidth  = 4;

smallSizeData   = 36;
mediumSizeData  = 48;
bigSizeData     = 60;

randomNumberStream = RandStream('mlfg6331_64', 'NormalTransform', 'Ziggurat');
subStreamNumber = 57162;
set(randomNumberStream, 'Substream', subStreamNumber);
RandStream.setGlobalStream(randomNumberStream);

%% Constants
UNIT8_MIN_VALUE = 0;
UNIT8_MAX_VALUE = 255;

%% Loading Data

mInputImage = imread('../RawData/Image0001.png');
mInputImage = mInputImage(:, :, 1); %<! Assuring Single Channel Image

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

%% Parameters
winRadius           = 1; %<! Loacl Window Radius
noiseProbability    = 0.05; %<! Probability of a Pixel to be affected by noise

%% Generating Noisy Image

numNoisyPixels  = round(noiseProbability * numPixels);

vNoisyPixelsIdx         = randperm(numPixels, numNoisyPixels);
vPepperNoisyPixelsIdx   = vNoisyPixelsIdx(1:floor(numNoisyPixels / 2));
vSaltNoisyPixelsIdx     = vNoisyPixelsIdx((floor(numNoisyPixels / 2) + 1):numNoisyPixels);

mNoisyImage                         = mInputImage;
mNoisyImage(vPepperNoisyPixelsIdx)  = UNIT8_MAX_VALUE;
mNoisyImage(vSaltNoisyPixelsIdx)    = UNIT8_MIN_VALUE;

%% Image Median Filtering

mFilteredImage = mNoisyImage;

for iCol = 1:numCols
    for jRow = 1:numRows
        winFirstRowIdx  = max(1, (jRow - winRadius));
        winLastRowIdx   = min(numRows, (jRow + winRadius));
        winFirstColIdx  = max(1, (iCol - winRadius));
        winLastColIdx   = min(numCols, (iCol + winRadius));

        vRowsIdx = [winFirstRowIdx:winLastRowIdx];
        vColsIdx = [winFirstColIdx:winLastColIdx];

        mLocalWin = mFilteredImage(vRowsIdx, vColsIdx);

        mFilteredImage(jRow, iCol) = median(mLocalWin(:));
    end
end

%% Display Results

hFigure = figure('Position', [100, 100, 900, 550], 'Units', 'pixels');
set(hFigure, 'Colormap', gray(256));
hAxes   = axes();
set(hAxes, 'Units', 'pixels')
set(hAxes, 'Position', [50, 50, numCols, numRows]);
hImageObject = image(mInputImage);
set(get(hAxes, 'Title'), 'String', ['Input Image'], ...
    'FontSize', titleFontSize);

hFigure = figure('Position', [100, 100, 900, 550], 'Units', 'pixels');
set(hFigure, 'Colormap', gray(256));
hAxes   = axes();
set(hAxes, 'Units', 'pixels')
set(hAxes, 'Position', [50, 50, numCols, numRows]);
hImageObject = image(mNoisyImage);
set(get(hAxes, 'Title'), 'String', ['Noisy Image'], ...
    'FontSize', titleFontSize);

hFigure = figure('Position', [100, 100, 900, 550], 'Units', 'pixels');
set(hFigure, 'Colormap', gray(256));
hAxes   = axes();
set(hAxes, 'Units', 'pixels')
set(hAxes, 'Position', [50, 50, numCols, numRows]);
hImageObject = image(mFilteredImage);
set(get(hAxes, 'Title'), 'String', ['Filtered Image'], ...
    'FontSize', titleFontSize);

%% Restore Defaults
set(0, 'DefaultFigureWindowStyle', 'normal');
set(0, 'DefaultAxesLooseInset', defaultLooseInset);

Enjoy...


FFT for long waveform I should take FFT of a very long vector of length 1,636,399,763 (i.e. length order is 10^9) preferably in MATLAB. However, on 128 GB RAM, my code throws Out of Memory error on line  part2Vec = repmat(part2Vec_OneRepetition, 1000, 1) from script below which  is attempting to build a vector of length 1,636,763,000 before even attempting to build totalVec which is of length 1,636,399,763 and take FFT which would be of length nFftLen = 2,147,483,648.
Any thoughts how do do this? For instance, can I take FFT in chunks and if yes, how I can do that?
%part1Vec = this is a "complex" column vector of length "1,636,763". Also, "complex" means each element is a complex number a+b*j.
%part2Vec_OneRepetition =  this is a complex column vector of length "1,636,763"
part2Vec = repmat(part2Vec_OneRepetition, 1000, 1); % This is a complex column vector of length "1,636,763,000".
totalVec = [part1Vec ; part2Vec]; % This is a complex column vector of length "1,638,399,763".
nLen = length(totalVec);% This is "1,638,399,763".
nFftLen = 2.^(nextpow2(nLen));
myFft = fftshift(fft(totalVec, nFftLen));

If compute time doesn't matter, you can overcome the memory problem easily with a dual for-loop implementation. But there's an FFT alternative with caveats: break up signal into M seqments, take M FFT's, add them up. The result is exact, limited in the lowest possible frequency it can represent. You could, however, do both; FFT for high frequencies, and for-loop for leftovers.

FFT in segments
Recall that FFT computes DFT, which in turn multiplies input by sines of cosines of varying frequencies. Let's take x of length 512, f=4 for left half, f=64 for right:

real(FFT(x)[1]), for example, was computed by multiplying:

And imag(FFT(x)[1]) would be computed by replacing the orange cosine of f=1 with orange sine of f=1. Now, what if N=512 is too long, and we decide to break it up into M=2 frames? Then, from those FFT's (FFT256) perspective, real(FFT(x)[1]) looks like

and from the FFT512's perspective, this looks like

which is real(FFT(x)[2]). In other words, FFT(x[:256])[1] + FFT(x[256:])[1] == FFT(x)[2].
Now imagine the highest frequency basis (sine or cosine) in either FFT256 or FFT512; what's the limit? It's [-1, 1, -1, 1], sample-to-sample (k=1/N). No matter how long the FFT is, this basis is always the greatest (normalized) frequency possible from x's point of view. As from above example, shortening segments corresponded to being unable to compute FFT(x)[1] - because that'd require a k=0.5 basis in FFT256:

It's important to understand "normalized frequency" to correctly combine the FFT's later. But basically, numerically the shorter segment computations are exactly the same as longer segment computations with higher-frequency bases.
For a billion point signal, the highest FFT basis frequency is k=500,000,000. What this actually means, is that the dynamic range of frequencies you can capture is 500,000,000:1 (look up normalized frequency). So, if in data you expect the ratio between the highest and lowest frequency to be 10 times less, then you can capture the entire spectrum by splitting up the billion points into 10 segments and adding up the results. If it's any more, you can capture a part of it with the above method, and another part with for-loops.

Additional notes:

If the input is real-valued, be sure to use rfft (real FFT), which takes half the memory.
Beware of any solutions proposing a window of any kind other than rectangular; they will distort the spectrum (but it may be practically insignificant).
Beware zero-padding
For very long DFT, be sure to use float64, and even that may not be enough depending on level of accuracy you seek

For-loops rDFT implementation: (Python but easy to port to MATLAB)
def dft(x):  # unnormalized
    N = len(x)
    reX = np.zeros(N // 2 + 1)
    imX = np.zeros(N // 2 + 1)

    for k in range(N // 2 + 1):
        for n in range(N):
            reX[k] += x[n] * np.cos(2 * np.pi * k * n / N)
            imX[k] += x[n] * np.sin(2 * np.pi * k * n / N)
    return reX - imX * 1j

Additionally
If working in Python, I recommend Numba for CPU-parallelism, and Cupy to run code on GPU. Both help with overcoming Python's slug for-loops and Numpy's single-threadedness.


Discrete signals in spectrum to Analog : DFT An analog signal $x(t)$ was sampled for a period $T$ of $1 \text{ms}$. The resulting, discrete signal is denoted as $x[n]$. Now you are given the spectrum $X[k]$ of the sampled signal $x[n]$.

Doubts regarding the attached images 1 & 2.

Please explain how this coefficients and signs for sin and cos term are considered??

For signal 1 ($x_1(t)$) :
The Fourier's transform ($F(f)$) of $Acos(2\pi ft)$ is $\frac{A}{2}(\delta(f-f_0)+\delta(f+f_0))$ and for $Asin(2\pi ft)$ it's $\frac{jA}{2}(\delta(f-f_0)-\delta(f+f_0))$

Replacing frequencies and amplitudes with the appropriate ones, the Fourier's transform of $x_1(t)$ is :
$$
F(f) = \delta(f-3000 Hz)+\delta(f+3000 Hz)+\frac{j}{2}\left(\delta(f-2000 Hz)-\delta(f+2000 Hz)\right)
$$

Hence the observed Diracs on fig. 1. You did not scale your axis according to the fft sampling while plotting it so the frequency does not appear clearly. The amplitudes are respectively 4 and 8 which is the right ratio so I guess you also did not scale the fft result (matlab doesnt do it itself), it should be 0.5 and 1 according to what I calculated.

For signal 2 it works the same way

edit : there is a mistake in the sampling frequency (1 kHz) you mention. To fit the graphs properly (and respect Shanon) it should be 8 KHz.


What's the use of the ripple in the Equiripple method and their effect in the filtering process? I read about the design methods of FIR filter which are:
windowing method, sampling frequency method and Equiripple method. And I don't understand the use of the ripple in the Equiripple method and their effect in the filtering process. Can anyone help me?

There is ripple to three different degrees depicted in this plot:

ripple is the variance that the filter gain is from the target gain.  in the above, there are 3 different frequency segments with 3 different target gains (well, two target gains are the same), and three different error weighting factors.

first frequency segement, from 0 to 0.2 Nyquist: target gain $-\infty$ dB gain (linear gain = 0, a "stopband"), about -67 dB maximum stopband gain (ripple = 0.00045).

second frequency segment,  from 0.2 to 0.6 Nyquist: target gain 0 dB (linear gain = 1, a "passband"), about $\pm \tfrac38$ dB deviation from 0 dB (ripple = 0.44).

third frequency segment, from 0.6 to 1.0 Nyquist: target gain $-\infty$ dB gain (another stopband), about -47 dB maximum stopband gain (ripple = 0.0045).

you will see that the ripple (in linear error terms) times the weighting (100, 1, 10) is about constant.

ripple is error.  it is the deviation of gain in the actual FIR filter from what your set target is.  this error with each band, besides being a tradeoff between each other (given a fixed number of FIR taps), it is a tradeoff with the steepness of the transition bands.  at least this is the case with the choice of FIR design method.  if the total "average" error is what is salient, then your design method might be least-squares (firls() in MATLAB) and if the maximum error (in each band) is what is important, then the design method is likely Parks-McClellan (firpm()).

now, theoretically, windowed FIR design using a good window (like Kaiser) can sorta compete with the two "optimized" and iterative methods above and there will also be a measure for ripple, but usually the ripple viewed as linear gain, is the same in the passband as in the stopband, but in dB is more in the stopband.  you can hardly see any passband ripple using a Kaiser-windowed design.  and the Kaiser-windowed designed FIR is maybe just a little longer than the optimal Parks-McClellan design, given the same transition band and stopband attenuation.

The "optimal" P-McC design will show ripple (whatever you allow in the tradeoff) in the passband and that might make the Kaiser-windowed design look a lot better to you.  Sometimes i pick a simple Kaiser-window design over P-McC, particularly if i want to upsample by just a factor of two.  In that case, if it's a windowed-sinc, i need not compute half of my output samples because i am copying them from the input (the P-McC does not have that feature).  but if i am upsampling by 4 or more, the P-McC will have shorter FIRs, enough so to make up for the requirement to compute the output sample instead of copying.


Hand tracking, state of the art I am looking for a framework for extracting hand region from the video feed from the webcam(analogous to Viola Jones HAAR detection of face).
This region of interest will serve as input for a CNN for hand gesture recognition, which I have already coded using python and CAFFE.
There are way too many frameworks and I am confused on which one to use,
priority would be to minimise the processing time. I already tried HAAR cascades for hand detection but they are pretty unreliable.
It would be most helpful if implementation is readily available in python.

Thanks in advance

The requirement for being able to detect the "profile" of hand gestures is not a straightforward one to satisfy. There are issues of form and rotation / scale independency which are not "easy".

There are two main problems here:

Segmenting the hands
Figuring out the gesture.

To segment the hands, let's assume that the hands are imaged over a constant single colour background. This means that the overall hand region can be isolated with a "simple" threshold.

The second problem is one of recognising form. In other words, what are the "shape" differences between "rock", "paper" and "scissors"?

Therefore, whatever is used should be "able" to discriminate between "holding two fingers extended" versus "all fingers retracted in the shape of a fist", etc.

The "easy" way out is to use template matching. Very "simply", take pictures of the hand gestures for each hand and for various orientations and use them to "search" for their existence in an image using cross correlation (or other similar concept in some feature space).

This however, done in the "simple" way is very costly (computationaly) and also scale dependent. So, if the hands are near the camera they are larger, if they are further away they are smaller and therefore very difficult to "catch" from static images that were acquired at some standard distance.

So, an alternative would be to use Morphological Operators and specifically Skeletonisation.

The main idea here is to apply skeletonisation on the extracted hand regions, quite a few times to end up with a set of "lines" describing the shape you are trying to recognise. In this way, in one pass you effect the skeletonisation and in another pass you collect the (very few, compared to the original) pixel locations that compose the gesture. To then figure out the gesture itself, you work with features extracted from the pixel locations.

The added benefit of this type of a solution is scale invariance. The relative positions of the "finger segments" of the "scissors" gesture are the same whether you hold it near the camera or far away from the camera (Two lines forming a V ending up in a big blob). With careful selection of the features and search process you can even use the same algorithm and training dataset for the left and right hands separately (for example, using graph search).

Hope this helps.


Why harmonic components appear only after a certain level when a signal is clipped? I recently observed this phenomenon that when a signal is clipped the harmonics start to appear only after a certain level.

The Python code to reproduce the effect is given below. The signal has 3 components at 30, 60 and 90 Hz. So a maximum possible peak to peak voltage will be 3. I clipped the signal at various levels and observed the magnitude spectrum but the harmonics start to visibly noticable only if the clipping goes below -2.5 and 2.5.

Is this a well-known phenomenon? What is the possible explanation behind this?.

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

labelsize   =   12
width       =   3.5
height      =   width / 1.618
lwidth      =   0.9

plt.rc('font', family='serif')
plt.rc('text', usetex=True)
plt.rc('xtick', labelsize=labelsize)
plt.rc('ytick', labelsize=labelsize)
plt.rc('axes', labelsize=labelsize)

def MagnitudeSpectrum(data, Fs, clip):
    P = 20*np.log10(np.abs(np.fft.rfft(data)))
    f = np.linspace(0, Fs/2, len(P))

    fig, ax = plt.subplots()
    fig.subplots_adjust(left=.20, bottom=.25, right=.96, top=.90)

    plt.plot(f, P, color='k', ls='solid', linewidth=lwidth, label='')
    plt.xlim(0, 500)
    plt.title('Signal Clipped BW '+str(-1*clip)+' and'+str(clip))
    plt.xlabel('Frequency [Hz]')
    plt.ylabel('Magnitude [dB]')
    # plt.show()
    fig.set_size_inches(width, height)
    fig.savefig('SignalMagSpectro_'+str(clip)+'.png', dpi = 600)

cliplevels = np.arange(2.2, 3.2, 0.1)

for clip in cliplevels:
    Fs = 44100.0
    t = np.arange(0, 10, 1/Fs)
    data = np.sin(2*np.pi*30*t)+np.sin(2*np.pi*60*t)+ np.sin(2*np.pi*90*t) +0.1*np.random.randn(len(t))
    data = np.clip(data, -1*clip, clip)
    MagnitudeSpectrum(data, Fs, clip)

Is this a well-known phenomenon?

Yes, of course. You will see harmonics as soon as your clip point is lower than the maximum amplitude in the time domain. The latter is a function of the relative phases between the harmonic components. In your case the max amplitude is indeed 2.5 (plus whatever the noise adds).

If you change the phases you will get a different clip point. For example, if you use cosine instead of sine, you will get a clip point of 3.0


How to convert multi-band Energy Time Curve (ETC) to Impulse Response (IR)? I am running my own sound simulator. Now I have obtained the Energy Time Curve (ETC, expressed by energy) by physical simulation, and they are stored in 8 discrete frequency bands (because decay of sound in the air is dependent on sound frequency). Now I am having problem getting the impulse response (IR, expressed by amplitude) of the measured system. I will need the IR to convolve with some audio files.

So I have $ETC(f,t)$, where $f\in\{f_1,f_2,...,f_8\},t=1,2,...N$. And the length $N$ is sufficient to capture the energy decay under a sample rate of 48kHz. How should I derive $IR(t)$ from the above $ETC$?

Thanks in advance.

You can't really get the impulse response from the ETC. This is because ETC represents the energy of the signal, while the impulse response contains time information too.
One may argue that ETC also has some time information but it's not the same. You need phase to be able to reconstruct the impulse response and this is simply not present in the ETC.
What you could possibly try to do is to sum the energy of the frequency bands to acquire the total energy (broadband), calculate the square root of the result to acquire the amplitude and add random phase, or even impose the phase characteristics you want. Keep in mind that this may lead to unstable results or complex values in your time series (IR).
In order to impose your desired phase response, you could convolve your signal (square root of total energy) with the IFFT of a "manufactured" signal that will have constant amplitude spectrum and the desired phase.
Once more, please keep in mind that arbitrary phase curves may lead to results that "do not make sense". No matter what, the resulting sound may, or may not sound natural/real. Most probably (haven't tested though) parts of it, such as the late reverberation tale may sound natural (it exhibits random-like phase characteristics) while the fist parts of the IR may not sound convincing.
The main issue with the manufactured phase is that real spaces exhibit mixed-phase characteristics. Some parts of the IR (early reflections) have minimum phase (most often) while late reverberation shows maximum phase characteristics.


Properties of convolution same as LTI systems? I have attached a snapshot from "signal processing first "
Where author says that properties of LTI systems are same as properties of convolution as shown in attached snapshot

Are these properties exactly same for both cases,convolution and LTI systems?and we don't need to apply any scaling etc?

The author simply means that the response of an LTI system $h[n]$ can be found by convolving the input sequence $x[n]$ with the $h[n]$. And in doing that all properties of convolution will hold true for LTI systems. Like, frequency response of an LTI system will simply become $X(e^{j\omega}).H(e^{j\omega})$, because convolution in time means product in $\omega$ domain.

You can look through all the properties on convolution operation at : ConvolutionProperties, and be sure that all of them will be held by an LTI System as well.

For example : Associativity : $$(x*h_1)*h_2 = x*(h_1*h_2)$$
In context with LTI systems, it means that if $x[n]$ is the input to two cascaded LTI systems $h_1[n]$ and $h_2[n]$, then we can actually, find out the equivalent LTI system by convolving the 2 cascaded systems as $(h_1*h_2)[n]$ and then find the output of this equivalent system by applying input $x[n]$ to it.

Similarly, you can see the relation between all of the convolution properties in context of an LTI system.


Repeated Fourier transform - what happens? I have a Fourier transformable complex function that is a function of independent real variable a.  Now I take the Fourier transform of it, giving me a complex function of real variable b.  Now I treat the resulting function as if it is in the original domain of a and again take Fourier transform of it - in stead of inverse Fourier transform as is usually done to get the original fucntion.  May be I will do this repeatedly.  What are relations between the original function and the transform function at various stages?

if you define the continuous Fourier Transform in a unitary manner, my preferred unitary definition is

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

$$ x(t) \triangleq \mathscr{F}^{-1}\{ X(f) \} = \int\limits_{-\infty}^{+\infty} X(f) \, e^{+j 2 \pi f t} \ df $$

then you can see a lotta symmetry and isomorphy in the forward and inverse transformation.  in fact they are exchangeable since $-j$ and $+j$ both have equal claim to squaring to be $-1$.  (i.e. if, in all of our textbooks and technical and math lit, every $j$ was replaced by $-j$ and vise versa, all of our theorems would be just as valid.  the choice of which imaginary unit to go with which direction of Fourier transformation is a convention.)

now, it's not hard to see where the duality theorem comes from.  Given the above, then

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

$$ X(-t) = \mathscr{F}^{-1}\{ x(f) \} = \int\limits_{-\infty}^{+\infty} x(f) \, e^{+j 2 \pi f t} \ df $$

from that, it's not hard to see that

$$ x(-t) = \mathscr{F} \Big\{ \mathscr{F}\{ x(t) \} \Big\} $$

and that

$$ x(t) = \mathscr{F} \Big\{ \mathscr{F}\{ x(-t) \} \Big\} $$

so it's not hard to see that if $x(t)$ is even symmetry ($x(-t) = x(t)$) then transforming it twice gets you back to the original.

and this

$$ x(t) = \mathscr{F} \bigg\{\mathscr{F} \Big\{\mathscr{F} \big\{ \mathscr{F}\{ x(t) \} \big\}\Big\}\bigg\} $$

so it's sorta like multiplying by $j$.  do it four times and you wind up with the thing that you started with.

you can use this fact to create an infinite number of Fourier transform pairs that are exactly equal to each other.  all you have to do is create the Fourier transform three levels deep and add each to the original.  each time you FT that, you get the same thing.

there are two simple functions that i can think of that have themselves as their own FT.  one is the Gaussian function and the other is the Dirac comb (both properly scaled).


Inverted values - digital resolution issue? first of all, I am not a technician, but I would really appreciate your helping with one (I hope simple or trivial) issue. I recorded a signal (skin conductance level - SCL, uS unit) with my device - 16 bit, range of the measured SCL is 1-20 uS.

In a recording software I changed digital upper range - "digital maximum" - from 1024 to 65536 (minimum = 0).

The issue is that when I record the signal and the signal reaches about the half of the physical maximum (the range the device can measure SCL = 1-20uS), from the recorded signal I see something like "inverted-across-zero" values, they have negative values but as it is from oscilatory parts obvious, when the signal is lower than the bound about 10,5 uS then it is above zero value, when it exceeds 10,5 uS it become inverted and the values are (deposited) under zero axis.

Can this "phenomenon" be caused due to my stupid change of digital maximum to 65536?
Thank you for any help!

M.

This looks like a "signed number" problem.

An unsigned binary number with n bits has a range from 0 to 2^n -1

But a Signed number weights the MSB as -2^(n-1) so the range is -2^(n-1) to 2^(n-1) -1

Consider an unsigned signal going from 0 to 255 being viewed as though it was a 2's complement number. Notice how it would look as though it flips from +127 to -128 and continue up to -1 (255 corresponds to -1).

Is that what you're seeing?

I notice that you're max digital value is 65536. But the maximum value a 16 bit unsigned number can hold is 65535 .. 2^16 -1


make a threshold frequency based on percentage of frequencies less than that? Using MATLAB signal processing toolbox, I have created signals like dirichlet, chirp and etc and and obtained their fourier and cosine fourier transform. Now I want to find the frequency where for example A% frequencies are less than that and use it as a thresholding frequency to compact the signal and eliminate unnecessary details?
Could you please tell me a little more about thresholding frequencies and how to choose them in signal processing?

A Matlab code and illustration is given below, with the  discrete cosine transform (DCT):

nSample = 512;
pDataThreshold = 90/100 ; % percentage, as a float in [0 1]

time = linspace(0,1,nSample)';
data = sin(2*pi*4*sqrt(time));
dataTrans = dct(data);
dataThreshold = quantile(abs(dataTrans),pDataThreshold);
dataTransThreshold = dataTrans;
dataTransThreshold(abs(dataTrans) < dataThreshold) = 0;
dataCompact = idct(dataTransThreshold);

plot(time,[data,dataCompact]);
xlabel('Time')
ylabel('Amplitude')
legend('Data',['DataCompact (',num2str(100*pDataThreshold),'%)'])
grid on


Distinguish harmonic sounds from non-harmonic I am trying to distinguish harmonic sounds in the range from $100\rm \ Hz$  to $1000\rm \ Hz$ in my audio recording.

I attach the spectrogram ($y$-axis is freqeuncy and $x$-axis is time) of a typical audio recording. The red square shows the harmonic sounds I want to distinguish, whereas the green one - the random sounds which I want to omit.

I can clearly see the harmonic components and consistent high-power fundamental frequency in the wanted sounds, while random sounds are not harmonic and frequency range is cluttered with low-power components.

My question is, how do I distinguish them?

I was thinking of calculating THD for both and selecting only the high results as harmonic sounds seem to carry more power, however, that's not always the case.

A standard difference between harmonic and non-harmonic sounds is that harmonic sounds tend to be periodic (in the time domain) in nature whereas non-harmonic sounds don't show any periodic properties since they're completely random.

One initial method you could try is simply performing autocorrelation in the time domain to see whether that portion of your signal is periodic. This could help you differentiate between the different sounds. There are definitely more complicated ways to test this. You could go to the extent of training a machine learning model based on certain features in the audio signal. But maybe you could just give the autocorrelation a try and see how well it performs for your application.