Comparing MFCC Features ,What do they represent? I know that MFCC features are the spectral envelope of the input signal but
I can't understand what do they mean and what do they represent  . and if I have two people saying the same word how can I compare the resulting features.

For example I have Mel Frequency Cepstral Coefficient  for the word "please":
first person:

second person:

I see these signals seem different but thay are representing the same word which is 'please' .

I've done these steps to extract the features:

framing with overlapping percentage equals to 50%.
Windowing : hamming window was used
perform  fft(signal);
take the squared magnitude of the ff(signal);
apply mel filter bank :number of filters = 26;
take the logarithm of  squared sum for each filter
perform dct .

Like already said, an eye inspection in this case is not possible since the data is just to high-dimensional. Thus what you want to do is to inspect the difference magnitudes of the frequencies per time frame.

I would suggest that you have a look at the cosine distance in this case. It's not the only way but a good possibility to compare MFCCs without distortion due to the different overall loudness of the word. The cosine distance is 1 - cosine similarity and is frequently used in text mining. In this area it might be important to be independent of the document length (document A with 100 words, 10 * "Hello" and 90 * "Bye" should be similar to document B with 200 words, 20 * "Hello" and 180 "Bye"). Another way of explaining this distance is perform a normalization step (project each data point on the unit sphere) and then measure the euclidian distance.

The cosine distance looks at the angle between two vectors (in your case you have one vector per time/frequency-bin) and therefore gets rid of the overall magnitude. This would by one way of comparing your MFCCs.


Can a Goertzel filter response be mapped to a Direct Form II Biquad IIR Filter I am looking at the structure of the Goertzel Filter, and it seems like I can map its coefficients to the standard implementation of the CMSIS DSP Biquad Cascade IIR Filter. Is this possible or should I give up and implement my own? The benefit here is the CMSIS Library can be had in a certified form for my application (I believe it just got rolled in). I believe I can compose the Goertzel response by cascading the Biquad with a Comb network, is this correct? I suppose my main point of confusion is centered from my inexperience in matters such as these, and if the Transposed Direct Form II of the cmsis biquad can demonstrate total equivalence to the Direct Form II suggested by this particular paper [edit]

try it, you'll like it.

the form in the Jocobson & Lyons paper is Direct Form II, with the coefficients explicit.  the form in the CMSIS source is a transpose of a Direct Form I, i believe, and the one in the Lyons & Bell paper is a Direct Form II (not transposed, but set $a_0=1$ because they should not have had that $a_0$ gain block in the path.).  the $b_0, b_1, b_2$ coefs come from the numerator and the $a_1, a_2$ come from the denominator and need to have their signs changed (since the scaled feedback signals are added and not subtracted in the feedback).

are you floating-point or fixed-point?  if the latter, the untransposed Direct Form I will serve you better.


How do I get a faster system response? I have this model in simulink (the graph is my output):

The step input has amplitude 0.5 m/s, and it steps up after 0.1 seconds. The gain $K_p=5$.
The saturation block is to keep the voltage between -3 and 3 volts.
The transfer function is the system, and it is given by $H(s) = \frac{16.94s + 579.5}{s^2 + 507.2s + 1224}$
The integration block is to convert velocity (in m/s) to position (in m). So the output is position.
I want my system response to be as fast as possible, and with no overshoot. On the graph above, an overshoot can be seen. To fix this, I added a $C_d=\frac{0.4s+1}{0.04s+1}$ term in the feedback loop.

As you can see there is no longer any overshoot. But it still takes about 3 seconds for the output to stabilize around the desired value.
How do I make the system response faster, such that the desired output value is reached sooner? Can I add a block in my simulink model perhaps?
Edit
Putting the $C_d$ block in the forward path sped up the response, but not by much. About 1 sec.

The transfer function is $H(s) = \frac{16.94s + 579.5}{s^2 + 507.2s + 1224}$
This transfer function has 2 poles, one slow pole at -2.4248 and a fast pole at -504.7752. The function has a slowish zero at -34.2. Good  news, your poles and zero are all in the left-half plane. It is much easier to control a system with zeroes and poles in the left-half plane than in the right-half plane. Your slow pole at -2.4248 will limit your performance. Your slowish zero at -34.2 will probably limit your performance by creating some overshoot, or maybe not. Your fast pole will have a limited impact on the closed-loop performance, so you can simply ignore it for now.
You could design your controller to cancel the slow pole at -2.4248 and the slowish zero at -34.2. This would be a lead lag compensator looking like this
$$H_{LeadLag} = K\frac{s + 2.4248}{s+34.2}$$
Then you could set a PI in series with the lead-lag compensator and tune the P and I gain to get the closed-loop performance you want.
$$H_{c} = H_{LeadLag} \frac{K_ps + K_i}{s}$$
That being said, I'm not sure my solution is particularly robust. You should check the Gang of four stability functions https://www.cds.caltech.edu/~murray/courses/cds101/fa02/caltech/astrom-ch5.pdf.
Edit 2 :
Another solution, simply use a PI controller and set the ratio of the P and I gain to 2.4248 to cancel your slow pole. This controller will not cancel your slowish zero, but maybe you don't need to cancel it to get the performance you want.
$$H_{PI} = K\frac{s+2.4248}{s}$$


How to find position, scale and rotation of a known shape/contour in an image? I have an image of an object.
The image shows a high contrast outline of the object.
It is guaranteed that the image shows exactly one object.
The shape/contour of that outline is known.
It is not simple (e.g. a circle) but not too complex either.
It is made up of (circular) arcs and straight lines.
Think of the outline of the head of mickey mouse in terms of complexity.
Given the image and the knowledge about the shape, how do I determine the position, rotation and scale of the shape and thus the object in the image?

I did some research and Hough transforms seem to be the way to go, but only for lines.
The wikipedia article states that

Altering the algorithm to detect circular shapes instead of lines is relatively straightforward.

I have the impression that as long as I can describe the shape mathematically, I can adopt the algorithm to that shape, because the mathematical model is what defines the transformation. Is that impression (while very rough) correct?
The article states further

For more complicated shapes in the plane (i.e., shapes that cannot be represented analytically in some 2D space), the Generalised Hough transform [12] is used, which allows a feature to vote for a particular position, orientation and/or scaling of the shape using a predefined look-up table.

Is my Mickey Mouse-shape a complicated shape?
Do I need the Generalised Hough transform?
Am I on the right track with this hough transform?

I found an implementation in openCV, which lacks scale detection:

finds arbitrary template in the grayscale image using Generalized Hough Transform
Detects position, translation and rotation

Basically speaking, I have an image and some shape description (think of a vector graphic of the outline) and I want to know where it is in the image. Does something like that exists?
I cannot share the image.

A Hough transform can certainly work. However, with position, scale, and orientation, you would have to use a 5-D Hough transform (2 DOF for position, 1 for scale, 1 for orientation), which might be computationally quite demanding.

As an alternative, you might want to look at SIFT or one of its extensions like SURF. It sounds like they can be adapted to do the job. And OpenCV has them.


What is the result of fourier transform of an image I am new to computer vision and Fourier transform concept. I have examined some codes. In those codes, an image is performed Fast Fourier transform by using the function np.fft.fft2(). As a result of invoking this function, we get a two-dimensional array filled with complex values.

What do these complex values refer to ?

Are they frequency values ?

A 2D FFT consists in comparing 2D data with a discrete set of evenly spaced 2D complex exponentials $e_{u,v}$ (made of a cosine for the real part, and and a sine for the imaginary part). Behind the scene, this comparison is exact: a comparison of a 2D data with enough 2D complex exponentials bears the same information as the original data, albeit easier to interpret.

For each complex exponential (as numerous as pixels in the image), this comparison (or correlation) yields one complex coefficient $c_{u,v}$. The modulus of $c_{u,v}$ tells you about the proportion of amplitudes related to complex exponential $e_{u,v}$; the phase provides information on how the complex exponential contribution is shifted in the image. All in all, this provides you with a novel representation of the 2D data, as a sum of shifted exponentials with different amplitudes.

So, the complex values are "complex amplitudes", corresponding to frequency values indexed by the indices $_{u,v}$.


Kalman Filter: Why do we decrease the state uncertainty regardless of the current measurement? I'm struggling with fully understanding the concept behind a Kalman filter.
For the sake of simplicity, let's ingore the input variable $u$ and assume constant process $Q$ and measurement noise $R$.
We can now predict the new state by using the state transition matrix $F$ and the old state estimation
$$x^-_k = Fx_{k-1}$$. Same does apply for the estimate uncertainty matrix : $P^-_k = FP_{k-1}F +Q$. Now we observe some data and depending on the ratio between the estimate uncertainty and the measurement uncertainty, we update our belief of the state estimation by calculating the innovation
$$x_k = x^-_k + K(z_k - Hx^-_k)$$
where $z_k$ is the current measurement, $H$ is the  observation matrix and $K$ is the kalman gain $P^-_kH^T(HP^-_kH^T+R)^{-1}$.
To incorporate the new measurement the state uncertainty is updated with $P_k = (1-KH)P^-_k$.

Do I see it correctly that $P_k$ decreases with each iteration because $ (1−KH)≤1$? Assuming that this is true, doesn't this mean that the uncertainty in the state estimation becomes very small with enough time steps and eventually tends towards $Q$ since the values within KH become smaller with each step?

If we are at a point where the values in $P$ are very small, which yes means that the uncertainty in our estimated state and the estimated measurement is very small, then the influence of the measurements is probably very small as the Kalman gain becomes small - I.e. even if the innovation $z_k-Hx^-_k$ is very large, we trust the prediction of the model, although perhaps the true observation would have been correct?

Suppose the values in $P_k$ are very small and we get measurements that are very different from the previous ones, but (for whatever reason) are true. The predictions of the kalman filter will deviate strongly from this and still have a small variance, simply because according to the previous data the prediction of the KL is more probable than the measurement? I.e. At a certain point  we have only predictions with a small variance, no matter whether they correspond to the truth and I see now way how the values in $P_k$ could grow again. But I think they should grow because if the innovation is big for some iterations, then there must be something wrong with the predictions

Thankful for any tip or reference

Do I see it correctly that Pk decreases with each iteration because (1−KH)≤1? Assuming that this is true, doesn't this mean that the uncertainty in the state estimation becomes very small with enough time steps and eventually tends towards Q since the values within KH become smaller with each step?

For constant $Q$ and $R$ it does decrease asymptotically, but it doesn't tend toward $Q$ because of the update step $P_k = (1-KH)P^-_k$.  Do a web search on "Steady State Kalman" for much more information than will fit into a StackExchange post.

If we are at a point where the values in P
are very small, which yes means that the uncertainty in our estimated state and the estimated measurement is very small, then the influence of the measurements is probably very small as the Kalman gain becomes small - I.e. even if the innovation zk−Hx−k
is very large, we trust the prediction of the model, although perhaps the true observation would have been correct?

If that is the case, then $Q$ or $R$ or both are wrong.  $Q$ and $R$ aren't chosen arbitrarily -- they must match the problem.  If (and only if) they do, then the Kalman is optimal in the least-squares sense.

Suppose the values in Pk
are very small and we get measurements that are very different from the previous ones, but (for whatever reason) are true. The predictions of the kalman filter will deviate strongly from this and still have a small variance, simply because according to the previous data the prediction of the KL is more probable than the measurement? I.e. At a certain point we have only predictions with a small variance, no matter whether they correspond to the truth and I see now way how the values in Pk could grow again. But I think they should grow because if the innovation is big for some iterations, then there must be something wrong with the predictions

Then $Q$ or $R$ or both are wrong, or the probability distributions aren't Gaussian and the assumptions of a plain old Kalman filter don't hold.
The Kalman filter assumes that the system whose states are being estimated is linear, that $Q$ and $R$ are known, and that the processes they describe are Gaussian.  If they aren't known (or are not matched to the problem at hand) then the Kalman won't give optimal results.  If the processes are not Gaussian then the Kalman filter will be the best linear filter in the mean-squared sense, but it won't be the best filter overall, and there may be some other sense (i.e., min-max) where it isn't optimal.
You seem to either be confused about the necessary condition (sometimes impossible to meet) that $Q$ and $R$ are known, or that the plain old Kalman assumes Gaussian processes.
If, for example, your measurement usually has a Gaussian error, but occasionally has a really bad measurement (i.e., it's got some almost-Gaussian long-tailed distribution), then you may be able to get close to an optimal filter with a filter that's a plain old Kalman at its core, but which monitors the error between $\hat z$ and the measured $z$, and if that error is much larger than predicted by $H\ P\ H^T$ you substitute in a much larger $R$ than "normal".
OTOH, if your system usually runs along smoothly, but every once in a while it acts like it got hit with a rock, you could again detect that $\hat z$ is not matching measurement, and substitute in a large $Q$ and recompute $P$ (especially if the "got hit by a rock" mechanism is known, and has a known effect on $P$).
Neither of these is indicative of problems with the plain old Kalman filter per se -- just that it is only optimal in one sense, and only for a certain model, and only when reality actually matches the model that you use to generate the filter.  Stray outside of the bounds where it is the right answer and -- you'll get the wrong answer.


Impulse Response Question For the differential equation
$$\frac{d^2y(t)}{dt^2}+\frac{3dy(t)}{dt} + 2y = x(t),$$

I was able to find the frequency response as $$H(j\omega)= \frac{1}{-\omega^2+3j\omega+2}.$$ However, I am not sure how to find the impulse response, which is the next part. Please help, and thank you in advance.

As Dilip pointed out in the comment above, you can get the impulse response using the inverse Fourier transform. However, a slightly easier method might be to use the Laplace domain instead; it's more amenable to easy inverse transforming via transform tables. First, recall that the frequency response is really just the $s$-plane transfer function evaluated along the imaginary axis. That is:

$$
H(j\omega) = H(s)|_{s=j\omega}
$$

Go through your frequency response expression and substitute $s$ back in for $j\omega$:

$$
\begin{align}
H(j\omega) &= \frac{1}{-\omega^2+3j\omega+2} \\
&= \frac{1}{(j\omega)^2+3j\omega+2} \\
\end{align}
$$

$$
H(s) = \frac{1}{s^2+3s+2}
$$

So now we have the transfer function in the $s$-plane. To get from here to an easily-inverse-transformable expression, we use a partial fraction decomposition on the transfer function (an online calculator is available here:

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

Given this form, we now refer to a table of common Laplace transforms to find:

$$
e^{-\alpha t} u(t) \Leftrightarrow \frac{1}{s+\alpha}
$$

where $u(t)$ is the Heaviside step function. The impulse response of the system is the inverse Laplace transform of its transfer function. Using the linearity of the Laplace transform and the above relation, we can therefore conclude that:

$$
h(t) = u(t) \left(e^{-t} + e^{-2t} \right)
$$


Do all systems in z domain are filters? Normally a system transfer function is represented by H(z)

I want to know wether all those system transfer functions are only representing a digital filter?

Or their any other thing/entity in z domain that can be represented by H(z)

One way to think about this:

Everything that can be represented as a rational function in the Z-domain can also be represented as a linear difference equation in the time domain. As such at can be interpreted (and implemented) as a filter. It may be non-causal or unstable, but it's still a filter.

You can certainly represent things in the z-Domain that are NOT rational functions, however, that would be a system transfer function anymore, because the transfer function is only defined for LTI (Linear Time Invariant) Systems


The frequency-domain spectrum obtained using MATLAB's FFT function seems inaccurate I am using MATLAB R2018b on a Windows 10 PC. In order to obtain the frequency domain information of a signal, I'm applying MATLAB's FFT function. The signal of interest is obtained from an experiment in which a Gaussian pulse is reflected by a  non-linear target. The reflected pulse is received on an antenna array. Thus, I expect the reflected signal's spectral component to show up on the frequency spectrum. However, the FFT produces what I can only describe as two impulses at the extremes of the frequency axis. I've attached the data being processed, and the results obtained. What exactly am I doing wrong? Are the results reasonable?

The input vector is a 100000x16 matrix (data from 16 different channels at 100000 samples each)

CODE:

% Code
Ns = 100000; %number of samples
f = linspace(0, 20e9, Ns); % frequency vector
frequency_spectrum = fft(time_domain_data);
plot(f, abs(frequency_spectrum)./max(abs(frequency_spectrum)));
title('Frequency Spectrum');
xlabel('Frequency (in Hz)');
ylabel('Normalized |Y(f)|');

The FFT produces values for both negative and positive frequencies.

Only plot the first half of the  FFT output.


Standard deviation of filtered noise In an other post i've found that the variance of filtered noise can be found by

$$\min_{\theta\in [\theta_1,\theta_2]}\left|H(e^{j\theta})\right|^2\sigma_x^2\le
\sigma_y^2\le\max_{\theta\in [\theta_1,\theta_2]}\left|H(e^{j\theta})\right|^2
\sigma_x^2\tag{1}$$

The derivation is based on spectral density functions and the relationship $$S_y(e^{j\theta})=S_x(e^{j\theta})\left|H(e^{j\theta})\right|^2\tag{2}$$ with $$\sigma_x^2=\frac{1}{2\pi}\int_{-\pi}^{\pi}S_x(e^{j\theta})d\theta=
\frac{1}{\pi}\int_{\theta_1}^{\theta_2}S_x(e^{j\theta})d\theta\tag{3}$$

The important point for me is that the variances $\sigma_x^2$ and $\sigma_y^2$ are related by a quadratic relationship.

I'm interested in the relation of the standard deviations $\sigma_y$ and $\sigma_x$, not the variances. So the question is: How is the standard devation $\sigma_y$ related to $\sigma_x$?

My naive idea (I'm yet not too experience with random variables in the frequency domain) is to take the square root and find  $\sigma_y= |H(\omega)|\sigma_x$. So the variances should have a linear relationship.

Is this correct or am I missing some point?

Thank's to MattL's comment I found the 'answer' myself... I was mixing up equations (1) and (2).

As MattL. pointed out, it makes sense to take the square root of equation (1), resulting in something like $$\min_{\theta\in [\theta_1,\theta_2]}\left|H(e^{j\theta})\right|\sigma_x\le
\sigma_y\le\max_{\theta\in [\theta_1,\theta_2]}\left|H(e^{j\theta})\right|
\sigma_x\tag{1}$$ if all terms are positive.

If one takes equation (2), applies the integral and the same scaling as in (3) the result is $$\frac{1}{2\pi}\int_{-\pi}^{\pi}S_y(e^{j\theta})d\theta=\frac{1}{2\pi}\int_{-\pi}^{\pi}S_x(e^{j\theta}) \left|H(e^{j\theta})\right|^2 d\theta \tag{2}$$

The left hand side equals the variance of y, so taking the square root gives the standard deviation. However the right hand side is different, as $\left|H(e^{j\theta})\right|$ is a function of frequency.

Hence in general  $\sigma_y\ne |H(e^{j\theta})|\sigma_x$.


how can we represent a Big Data set in Graph Signals. (Big Data in Graph signal Processing) Hi

I have read this article https://pdfs.semanticscholar.org/65d6/1afd9c35b0a75d9de77c2a4a2428af0f7f7b.pdf about Big Data analysis in Graph Signal Processing

I have a couple of questions :

first of all, consider a big dataset which is given to us, it has 150 rows and 365 columns which are sensor measurements and days of one year respectively($D$).

In Example Application section it tries to represent the big dataset in the product of two graphs (the article introduces 3 types of products but I think the dataset have been tested on the Cartesian product).
It breake the data set into two graphs (time series($T$) with 365 nodes and the sensor network($A$) with 150 nodes).

Because there was no tag named Graph signal processing, I could not mention it.

the article has represented this two graph in FIG1 a,b.

Now my question :

if I correctly understand!

if there was no time we just need to decompose the sensor matrix ($A$) and after that, we could take Fourier transform, now during the time we just have more samples from our sensors so why we should eig-decompose our product($T \cdot A$) instead of sensor network($A$).
I mean how increasing the number of sensor data from our sensors, represents itself in graph product.

or why matrix $D$ can be broken to ${T, A}$?(and definitely ${T,A}$ can be processed parallelly)

I would appreciate your solutions.

Just read it myself.  The justification for using the term, "Fourier Transform" comes from the idea that for a discrete time periodic time series, convolution reduces to a sequence of circulant shifts and inner products.  The "filters' correspond to diagonal matrices, so when he forms his data structures, the matrix eigendecomposition  is "equivalent" to a DFT.  The circulant shift matrix is the justification for using the word "Fourier".  When I was reading through the paper, I recalled that there is something called spectral graph theory where you do SVD (or orthogonal decompositions) on adjacency matrices, and the singular values (or eigenvalues) have some sort of significance.  This paper appears to show how graph products where one of the graphs corresponds to periodic sequences has meaningful Fourier interpretation. There are papers on graph polynomials as well.

The example didn't particularly enthuse me with a burning desire to master the subject, but I have notices a growing literature associated with anomaly detection. I have a tendency toward cynicism which is often misguided.

To answer the question that I think you asked:

You can do an orthogonal decomposition on any adjacency matrix ( directed, undirected, acyclic, and with cycles) but to be a proper DFT abstraction, you need a cyclic matrix (graph) someplace.

The only concern I had with the paper (other than my partial understanding) was with the efficiency claims associated with sparse matrices.  Sparse matrices involve a lot pointer chasing and modern parallel vector processors that do dense linear algebra very well are not that great with sparse linear algebra.  Thread based architectures like GPU's seem better.


How to stitch together MFCCs from multiple frames? I am trying to record a sound two times and compare the two recordings. I am playing exactly the same sound both times (for example, play a small .wav file twice and record it each time), so I hoped that my attempt would not rise to the complexity of Speech Recognition...

I have tried computing the DFT (using FFT) and comparing the two power spectrums of the signals (more precisely, comparing their MSE - Mean Squared Error), but the results were discouraging (sometimes for recordings of the same sound the error was greater than for different sounds)...

Then I found out about the MFCCs and computed them for each signal. That leaves me with 12 coefficients for each frame I devided my signal into (details below). How should I combine these groups of 12 coefficients to get some value(s) that describe the whole signal?

(Appending them and extracting the MSE for the two arrays did not return encouraging results either...)

Below are some details about the way I computed the MFCCs and values I chose for computing it. Maybe it's here that I'm doing something wrong:

(sample rate = 22000 Hz)

Pre-emphasis.
Framing with frame size = 512, frame overlap = 200. (because "each frame should be ~20-30ms long" and 512 / 22000 = 0,023 seconds = 23 ms. As "the overlap should be of ~10ms and 200 / 22000 = 0,009s = 9ms)
Apply Hamming window to each 512-size frame.
Apply FFT to each windowed 512-size frame => 512 magnitudes for my signal, from witch I use only the first 256. Domain: Hz from 0 to 22000 (to 11000, respectively)
Compute the Mel Filter Bank: (min frequency = 300 Hz, max frequency = 11000 Hz)

Compute mels from min and max frequences and then compute 26 equally distanced values between these two mel values.
Convert them back to Hz => array of 28 frequency filters.
Compute a filterbank (filter triangle) for each three consecutive values => 28 - 2 = 26 filterbanks.
Pass the whole power spectrum (256 magnitudes) from step 4 through each triangular filter to get a "filterbank energy" for each filter => 26-sized array of energies.

Apply log: ln(each energy) (they are still frequency-domain values)
Apply DCT to the logged energies => 26-sized array of (time-domain?) values.
Take only the first 12 => 12 MFCCs for each 215-sized frame of my signal.

How should I combine these groups of 12 coefficients to get some
  value(s) that describe the whole signal?

I don't think "flattening" the sequence of $215 \times 12$ coefficients into a $1 \times N$ vector, and then computing the distance between these vectors to compare two sounds is a good idea. This approach collapses the temporal information - the way a sound evolves in time. For example, let us say that you decide to use the standard deviation and mean of each coefficient to "summarize" the sequence of MFCC vectors (this is a technique used in some applications). If you use this approach, both a sound and the same sound played backwards will yield exactly the same feature vector! This might work for simple tasks like recognizing different music instruments, but as soon as the temporal information needs to be preserved (distinguishing "halo" and "hola", or a drum instrument playing "boom boom tchik" vs "boom tchik boom"), thinking in terms of vectors rather than sequences no longer works.

This might not be the answer you expect, but I suggest you to stop reasoning in terms of vector comparison, and rather think of a way of comparing the sequences of MFCC frames themselves. Your approach - which is akin to computing the raw pairwise distance between sequence of vectors - is not robust to difference in timing between the two sequences. The solution is to use Dynamic Time Warping. It attempts to optimally align the two sequences in time (so that comparing "hoooolaa   " with "   holaaa" yields the same score as comparing "hola" and "hola"). This technique is good enough for basic voice command recognition.

Another approach would be to model each sound by a HMM with a few states; and to use the distance between models as a similarity metric. This can be more computationally efficient, and put you in control of how strict you want to respect temporal information (by increasing the number of states).


Histogram matching of two images using cdf I have already calculated the cdf values of each image now i have to map the reference image values to target .
I cannot use histogram equalization or interpolation
I simply have to map them .

function S = histM(t, target)
 if size(target,2)>1
   %an image, replace the image with its cdf
   target = getImageCDF(target);
 end

%get the CDF for the input image
 tcdf = getImageCDF(t);

% .... what shall be done here to map the values of reference to target

%subfunction: compute the CDF for an image that started as uint8
 function cdf = getImageCDF(img)
 bins = 0:255;
 H = hist(img(:), bins);
 Hmod = H + eps(sum(H));
 cdf = [0, cumsum(Hmod)/sum(Hmod)];

First you'll have to compute the histogram of one of the two images.

H = hist(img(:), bins);

Next find the cdf of the image:

 cdf = [0, cumsum(Hmod)/sum(Hmod)];

Next you'll have to make the second image follow the exact same cdf of the first image. This is usually known as histogram specification.
Here is the algorithm:

Say you have a 3 bit image, that is $8$ different intensity values $0 > 7$.
Now as for the second image

Find it's cdf as well.

Now as seen in the following image. The $r_k$ is the pixel value of the second image.

$s_k$ is the cdf of the second image.

$z_k$ is the cdf of the firsk image (desired histogram).

$z_{k}*$ and $n$ are computed as the following:

$(z_k - s_k) >= 0$ for smallest value of $n \dots$ and $n = g_k$

How exactly:
Let's compute $z_k*$

$z_k - s_k = 0 - 0.15 = -0.15 >! 0 $ (not greater than zero)
increase $z_k$ and keep $s_k$ as is until you satisify the criteria.

Until you reach:
$0.15 - 0.15 = 0 >= 0$ (check)
Then $z_k* = z_k$
$n = r_k $

Here are some illustration images:


Downsides of higher order filtfilt's I understand that with a standard filter, higher order will yield sharper frequency selection at the cost of increased phase delay. However, with so-called zero-phase-distortion filters like that provided by scipy.signal.filtfilt, will a higher order get you sharper frequency selection "for free"? (I know higher order filters take more computation time, but let's ignore that)

The sharper the filter is in the frequency domain, the longer the impulse response will be. This typically leads to "time blur" or "ringing" in the time domain. In addition, a zero phase filter is non-causal, so you get "pre-ringing" and any sharp onsets or transients in the time domain get degraded.

The long impulse response also leads to a long "transition" region at the beginning and the end of the signal. As a very rough rule: if you have a filter that does something drastic over 1Hz of bandwidth, it will need about 1 second of time domain signal until the effect of the filter is fully established.

The optimal choice is very dependent on the specific requirements of your application.


Licence plate enhancement Is it possible to enhance the image below and make it more readable?

I have tried apply binarization but there is some kind of shadow.

So I have tried to convert image to gray scale and then equalization

It is better, but still there is a problem with shadow and international vehicle registration code is not readable.

I guess image denoising and image super-resolution techniques could be used to enhance license plate images.
Here are a few relevant links:

Image denoising based - SNIDER: Single Noisy Image Denoising and Rectification for Improving License Plate Recognition
Image super-resolution based - Multiframe Superresolution of Vehicle License Plates Based on Distribution Estimation Approach


Compensating the Group Delay of an Analog Filter using DSP I am having a problem in my current project where I need to measure different biological signals (like ECG) simultaneously.

As the instrumentation hardware is developed by different people over time, they are not identical to each other. All of them has different orders of low-pass and high-pass filters before their final amplification stage. These filters are mostly concentrated around the passband of [0.5-35] Hz. I want to compensate for the phase dispersion caused by those analog filters by designing my own digital filters. I have measured the transfer functions of those channels using chirp signals, to get an idea, please check the following figure. There might be small numerical errors in the calculation of the transfer function.

Not only eliminating the phase dispersion in the individual channels but also compensating the different amount of delay of those channels is important so that I can analyze fully "synchronized" signals. So far, I have been trying to design an all-pass filter in Matlab using iirgrpdelay function. However, the created compensation filters are not satisfying. The reason for that might be the huge variation in group delay (see the first figure) between [1 15] Hz, the region I want to compensate. Increasing the order of the filter doesn't help since I mostly get the following error from Matlab:

Poorly conditioned Hessian matrix. Cannot accurately compute the optimization because either the approximation error is extremely small (try reducing the
number of poles or zeros) or the filter specifications yield huge magnitude variations, such as mag=[1 1e9 0 0].

I don't want to leave the question this broad but I am asking for suggestions to design such compensation filters in the digital domain so that I can neutralize this phase dispersion and delay in the analog instrumentation.

It's often easier to design FIR filters for compensating group delay. At the same time they could also compensate the magnitude if necessary. The easiest method is to use a complex least squares method, which boils down to solving a system of linear equations for the filter coefficients.

The difficult part is to choose an appropriate bulk delay in the definition of the desired phase response. A good starting point is to choose the average delay to be close to half the filter length.

Take a look at this answer for an example of FIR group delay equalization.

You can use the Matlab/Octave function lslevin.m for designing FIR filters with prescribed magnitude and phase responses.


Are there any cases where the negative frequencies of an FFT are used? In performing FFTs and looking at spectra for real signals (say, taken from a sensor) we ignore the negative frequencies. I understand that the negative frequencies often reflect the positive frequencies because their imaginary components are complex conjugates. Is this only the case for real signals? Would we be able to notice a difference if imaginary components of a signal did not cancel out (as in hearing an audio signal)?

A strictly real vector in one domain implies a complex conjugate mirrored vector in the other domain, and vice versa, for any FFT.

The "negative" frequency results of an FFT are very commonly used in SDR (Software Defined Radio), where the signal of interest might consist of only negative frequencies (LSB or lower sideband SSB is one obvious case) in a baseband IQ stream (after quadrature downconversion from strictly real higher frequency RF).  Or the SDR software may display a wideband spectrum showing signals both above and below the IQ baseband center frequency, which are displayed from an FFTs negative frequency bins thru positive frequency bins.

Some USB SDR devices (Funcube, et.al.) communicate the complex IQ data stream as stereo audio inputs to the computer, so one might actually be able to listen to the complex IQ data as a stereo left and right pair of audio (after appropriate sample rate conversion, if needed).


Calculate average from a power spectrum I want to calculate the average of a power spectrum in the complete frequency range. I have done a stft and want to calculate power in dB and then the average. Please let me know from the code below if my understanding is correct and how to calculate the average of the power spectrum.

y, sr = librosa.load('Noise Recorder.wav',sr=None,duration=2)
n_fft = 8192
S = librosa.stft(y=y,n_fft=n_fft,win_length=n_fft,hop_length=n_fft,window='hann')     # Short FFT
F = librosa.fft_frequencies(sr=sr, n_fft=n_fft) # Frequency components
Abs = np.abs(S)  # Same as mag, phase = librosa.magphase(S)
Abs_avg = np.mean(Abs, axis=1)
Power_avg =librosa.power_to_db(Abs_mean**2) # dB
Scale_dB = 86.85    # Scaling factor

Power_Scaled_dB = Power_avg + Scale_dB  #Plot #Max:114dB at 1KHz
RMS = np.sqrt(np.sum(Power_Scaled_dB**2) / n_fft) # n_fft=8192 #Ans:24.4

###### Additional RMS tryouts #####
# Method 1 : Based on PSD formula
Px = 10**(np.max(114/10)) # Ans: 251188643150.9
RMS_Px = np.sqrt(np.sum(Px**2) / n_fft) # Ans: 2775268639.5
RMS_Px_dB = 10*np.log10(RMS_Px) # Ans: 94

# Method 2: Based on amplitude fomrula
A = np.sqrt(Px) # Ans: 501187.2
RMS_A = A / np.sqrt(n_fft) # nfft= 8192(fft length), Ans: 5537.388930343294
RMS_A_dB = 20*np.log10(RMS_A) # Ans: 74

With a 114dB 1Hz Sine wave input signal, after stft i obtain the psd in the plot
Now I want to derive the PSD from the spectrum, to further calculate the average/RMS of the spectrum. Given by the relation PSD and from there I want to calculate RMS given by the relation RMS/Leq. I am unsure how to derive the PSD from the Power_avg
 and

By Welch's method, you can calculate the power spectrum by averaging the magnitude of a bunch of FFT frames (or time frames of the STFT). If your signal is $x(t)$, and its STFT is $X(\omega, \tau)$, where $\omega$ is the frequency bin, $\tau$ is the time frame and $T$ is the total number of frames, the PSD is
$$P_x(\omega) = \frac{1}{T} \sum_{\tau=0}^{T-1} |X(\omega, \tau)|^2 $$
When calculating the STFT, you need to make sure you get constant overlap add. This means, with the Hann window, you need a hopSize of N/2 (where N is your window length). I'd also keep the window length shorter, about 256/512 samples. This is not compulsory in your case, because your signal is stationary. Had it not been stationary you should have chosen a smaller window length. Here is your code, modified.
    y, sr = librosa.load('Noise Recorder.wav',sr=None,duration=2)
    n_fft = 8192
    win_len = 512
    S = librosa.stft(y=y,n_fft=n_fft,win_length=win_len,hop_length=win_len/2,window='hann')     # Short FFT
    F = librosa.fft_frequencies(sr=sr, n_fft=n_fft) # Frequency components
    Abs = np.abs(S)  # Same as mag, phase = librosa.magphase(S)

    PSD = np.mean(Abs**2, axis=1).  #this is the power spectrum
    PSD_dB = 10*np.log10(PSD). #power spectrum in dB
    RMS = np.sqrt(np.mean(PSD**2))) #RMS of the power spectrum
    RMS_dB = 20*np.log10(RMS) # RMS of power spectrum in dB

    ##Theoretical RMS (in dB) for a sine wave of amplitude A dB
    A_lin = 10**(A/20.0)           #convert dB to linear amplitude, A = 114 in your case
    RMS_exp = 20*np.log10(A_lin/n_fft)   #expected RMS in dB

`


Adding encryption to a ham radio Suppose I want to encrypt the voice data before sending using a ham radio, for example this unit.

Assume the encryption is simply AES encryption/decryption using predefined fixed key.

I am thinking of the following implementation:

1) Between the radio's transmitter/receiver and its logic unit, I'll place my circuit.
2) Low pass filter the signal then ADC it (with 10-bit accuracy, but this is irrelevant) with sampling frequency identical to that used by the radio. Ideally, this rate will already be higher than the max. freq. so as to avoid aliasing.
3) Apply the AES algorithm to each 128-bit segment of data.
4) DAC.

Going forward, I would implement some digital signature methods, thus signing the hash of the data, every 100 128-bit segment (or any other number). This will add 512-bit of data every 100*128 bits.

My question: 1) Is there any comment on the implementation, from a signal processing / communication point of view? 2) The added 512-bit signature will decrease the bit rate by (12800-512)/12800 => 4.16%. Is this a problem?

Digital audio data that is transmitted wirelessly is usually encoded in some efficient format. It sounds like you're talking about broadcasting raw PCM samples that are encrypted before transmission. That's certainly unusual.

Unless you are already following some kind of standard that I'm not aware of, you're probably best suited by packing your audio data into the payload some higher-level protocol that includes synchronization data to recover framing, the encryption that you desire, and most likely error correction. Your approach of just encrypting raw PCM samples with no error correction has the unpleasant side effect that a single bit error in the encrypted data will turn your entire block of samples into garbage.


LMS algorithm for modeling step-size ambiguity! Behrouz Farhang-Boroujeny in his adoptive filters, 2nd ed., p. 155 ,told:

  It is sufficient for stability: 1/(3*tr(R)).

But his book have attached mfile for modeling:

%   Modeling
% Last updated on April 28, 1998
%
itn=input('\n No. of iterations?      ');
sigman2=input('\n Variance of the plant noise?      ');
sigman=sqrt(sigman2);
wo=input('\n Plant impulse response (vector, w_o)?      ');
a=size(wo);
if a(1)<a(2)
   wo=wo';
end

N=input('\n Length of the model (N)?      ');

h=input('\n Coloring filter impulse response (vector, h)?     ');
a=size(h);
if a(1)<a(2)
   h=h';
end

Misad=input('\n Misadjustment (e.g., 0.1 for 10%) ?     ');
mu=Misad/(N*(h'*h));

You can see:

mu=Misad/(N*(h'*h))

And it is not matching to his formula in his book ...
I can't understand where the hell it is coming from?

More information about problem:

Scheme of problem:(Note the source of image is irrelevant to my question, it's only image)

Core of algorithm: (Note, it is average over runs for smooth learning curve)

for k=1:runs
    x=filter(h,1,randn(itn,1));% the randn creates input sequence that colorize by filtering through filter with h impulse response
    d=filter(wo,1,x)+sigman*randn(itn,1); % Create desired by sum of noise and y(main plant output)
    w=zeros(N,1);

    for n=N:itn
        xtdl=x(n:-1:n-N+1);
        e=d(n)-w'*xtdl;
        w=w+2*mu*e*xtdl;
        xi(n)=xi(n)+e^2; % Creation of learning curve, gradient(e^2)=-2e(n)x(n)
    end
end
xi=xi/runs;

Note: 1/(3*tr(R)) make algorithm unstable ...!!! But without colorizing filter this will converge.

I have the first edition of Behrouz Farhang-Boroujeny's Adaptive Filters book.  I found it useful and it was definitely more practical in terms of implementing adaptive filters than other textbooks like those from Haykin and Sayed, primarily because of the included Matlab code.  However, like any topic in the area of adaptive filtering, I would use it with a grain of salt.

There are numerous variations of the LMS algorithm, primarily because it is readily implemented and not computationally expensive compared to other adaptive algorithms such as the Recursive Least Squares (RLS) algorithm.

I think a lot of the LMS variations you can find (particularly in Farhang-Boroujeny's book) come about because people tweak one of the more traditional variations, ending up with something that performs well for their particular problem.  These often make it into journals because they are different and are indeed novel implementations.

I would suggest you read about the standard LMS algorithm and the normalized version of this algorithm, the NLMS.  Once you feel comfortable with your understanding of these algorithms, then implement them and apply them to your problem.  Tuning the step-size can then be done by trial-and-error.  If you are keen to understand the mathematical details of how the step-size is determined, I suggest the Haykin and Widrow textbook.  It has in depth analysis of the convergence behavior of LMS-based algorithms.

Here is an example (my code) of the LMS algorithm in Matlab.

function [ prediction_error, weights ] = LMS_Algorithm( regressive_sequence, ...
    step_size, number_of_taps )

% This script-file implements the Least Mean-squares (LMS) adaptive
% algorithm.
%
%
%% See Also
% LMS_Algorithm, NLMS_Algorithm, RLS_Algorithm, APA_Algorithm

% Author:  Michael R. Wirtzfeld
% Modification Date:  Monday, February 27, 2012
% Creation Date:  Wednesday, February 15, 2012

weights = zeros( number_of_taps, size(regressive_sequence, 2) );
prediction_error = zeros( size(regressive_sequence) );

for index = (number_of_taps + 1):1:numel( regressive_sequence )

    tap_inputs = [ regressive_sequence( (index - 1):-1:(index - ...
        number_of_taps) ) ]';

    prediction = weights(:, index)' * tap_inputs;

    prediction_error(index) = regressive_sequence(index) - prediction;

    weights( :, (index+1) ) = weights( :, index ) + ...
        step_size * tap_inputs * prediction_error(index);

end;  % End:  for index = 1:1:numel( regressive_sequence )

weights = weights( :, 1:1:numel( regressive_sequence ) );

%% References

% Section 9.5, Summary of the LMS Algorithm, Adaptive Filter Theory, Third
% Edition, S. Haykin, Prentice Hall, 1996

Here is an example of the NLMS algorithm in Matlab.

function [ prediction_error, weights ] = NLMS_Algorithm( regressive_sequence, ...
    step_size, number_of_taps )

% This script-file implements the Normalized Least Mean-squares (NLMS)
% adaptive algorithm.
%
%
%% See Also
% LMS_Algorithm, NLMS_Algorithm, RLS_Algorithm, APA_Algorithm

% Author:  Michael R. Wirtzfeld
% Modification Date:  Monday, February 27, 2012
% Creation Date:  Wednesday, February 15, 2012

a = 0.0001;

weights = zeros( number_of_taps, size(regressive_sequence, 2) );
prediction_error = zeros( size(regressive_sequence) );

for index = (number_of_taps + 1):1:numel( regressive_sequence )

    tap_inputs = [ regressive_sequence( (index - 1):-1:(index - ...
        number_of_taps) ) ]';

    prediction = weights( :, index )' * tap_inputs;

    prediction_error(index) = regressive_sequence (index) - prediction;

    normalization_term = a + norm(tap_inputs)^2;

    weights( :, (index+1) ) = weights( :, index ) + ...
        (step_size / normalization_term) * tap_inputs * prediction_error(index);

end;  % End:  for index = 1:1:numel( regressive_sequence )

weights = weights( :, 1:1:numel( regressive_sequence ) );

%% References

% Table 9.2, Summary of the NLMS Algorithm, Adaptive Filter Theory, Third
% Edition, S. Haykin, Prentice Hall, 1996

These functions were used as part of an investigation to see how a standard tap-delay line implementation of the LMS and NLMS algorithms could model a second-order auto-regressive sequence.

As Fat32 said, most of the textbooks on the topic of adaptive filters use different notations and suggest different recommendations regarding how to set parameters, so it's no wonder that you feel uncertain about what to do.

I hope this helps.


Does the ROC (Z-Transform) always start/end from/in a pole? Can I assume that the R.O.C. from a Z-Transform will always start from or end in a pole or 0/infinity?

Using an example:

$$H(z)=\frac{(z+1)(z-1)(z+j)(z-j)}{z^4}\space,\space\space j=\sqrt{(-1)}$$

So we have 4 zeros in the unit circle equally spaced by a $\frac{\pi}{2}$ angle and 4 poles in the origin ($z=0$). If not specified, can I assume the the ROC wil necessarily start from exclusive zero and go towards infinity? Or maybe possible to exist a ROC which is a ring in the Z-plane (which will not start/end in a pole)?

Since the value of Z-transform tends to infinity at poles, poles can not be a part of ROC. So ROC either starts from a pole or end in a pole. For example,

If $x[n]$ is a right-sided sequence, then the ROC extends outward from the outermost pole in $X(z)$
If $x[n]$ is a left-sided sequence, then the ROC extends inward from the innermost pole in $X(z)$
If $x[n]$ is a two-sided sequence, the ROC will be a ring in the z-plane that is bounded on the interior and exterior by a pole
If $x[n]$ is a finite-duration sequence, then the ROC is the entire z-plane, except possibly $z=0$ or $|z|=\infty$

From all these cases, one can say that poles decides the boundary of ROC.

The example you provided is a finite-duration sequence and its ROC is entire z-plane except at z =0.


Periodicity of the discrete-time Fourier Transform The DTFT of a sequence $x[n]$ can be written as $$X(e^{j\omega}) = \sum_{n = -\infty}^{\infty} x[n] e^{-j\omega n}.$$
Is the smallest (fundamental) period in frequency of the DTFT always $2\pi$? Or can it be smaller than $2\pi$?

In addition, I was wondering why we cannot use that same notation for the argument with the CTFT. We typically denote the continuous-time Fourier Transform by $X(j\omega)$. But what's wrong with using $X(e^{j\omega})$ where $$X(e^{j\omega}) = \int_{-\infty}^{\infty} x(t) e^{-j\omega t} dt?$$

I know for a fact that this is wrong because the notation inherently implies the periodicity of $X(e^{j\omega})$, since $X(e^{j(\omega + 2\pi)}) = X(e^{j\omega})$, which is obbiously not true in this case. But I don't understand what mathematically does not permit us to use $e^{j\omega}$ as an argument in the continuous-time case.

The DTFT is always $2\pi$-periodic. However, it can also have a smaller period, namely a fraction of $2\pi$. Take any sequence $x[n]$ for which the DTFT exists and insert $L-1$ zeros between the samples. The DTFT of the new sequency $\hat{x}[n]$ can then be written as

$$\hat{X}(e^{j\omega})=\sum_{n=-\infty}^{\infty}\hat{x}[n]e^{-jn\omega}=\sum_{n=-\infty}^{\infty}\hat{x}[nL]e^{-jnL\omega}\tag{1}$$

$\hat{X}(e^{j\omega})$ as given by $(1)$ clearly has a period of $2\pi/L$.

Concerning your second question, the reason why a continuous-time Fourier transform $X(j\omega)$ cannot be written as $X(e^{j\omega})$ is that

$$X(j\omega)=\int_{-\infty}^{\infty}x(t)e^{-j\omega t}dt\neq \int_{-\infty}^{\infty}x(t)(e^{j\omega})^{-t}dt\tag{2}$$

simply because generally

$$e^{-j\omega t}\neq (e^{j\omega})^{-t}\tag{3}$$

unless $t$ is an integer (as is the case with the DTFT).

Note that if $(3)$ were not true,  i.e., if $e^{-j\omega t}= (e^{j\omega})^{-t}$ were true, we could easily show that $e^{-j\omega t}=1$ for all $t$. Just write $\omega=2\pi/T$ and you get

$$e^{-j\omega t}=e^{-j2\pi t/T}=\left(e^{-j2\pi}\right)^{t/T}=1^{t/T}=1$$

which is of course absurd.

Proof of Eq. $(3)$:

Let $z$ be a complex number with $|z|=1$ (the magnitude is irrelevant here):

$$z=e^{j\theta}\tag{4}$$

Let $a$ and $b$ be real numbers. Then

$$z^{ab}=e^{j\theta a b}\tag{5}$$

And with $z^a=u$

$$\left(z^a\right)^b=u^b=e^{j\arg\{u\}b}\tag{6}$$

With $\arg\{u\}=\text{pv}\{\theta a\}\in (-\pi,\pi]$, where $\text{pv}$ denotes the principal value, it is clear that generally

$$z^{ab}=e^{j\theta a b}\neq \left(z^a\right)^b=e^{j\,\text{pv}\{\theta a\}b}\tag{7}$$

There are two cases where $(5)$ and $(6)$ are equal:

if $\text{pv}\{\theta a\}=\theta a$, which is the case if $\theta a\in (-\pi,\pi]$.
if $b$ is integer, since

$$\text{pv}\{\theta a\}=\theta a+2\pi k\tag{8}$$

with some appropriately chosen integer $k$ (such that the result is in the interval $(-\pi,\pi]$), then

$$\left(z^a\right)^b=e^{j\,\text{pv}\{\theta a\}b}=e^{j(\theta a+2\pi k)b}=e^{j\theta ab}=z^{ab},\qquad b\in\mathbb{Z}\tag{9}$$

[Note that $(9)$ would also hold for rational $b$ as long as $kb$ is integer.]


Do DCT results contain phase spectrum? I've been working much with numpy.fft lately where just like the documentation says:

  When the input a is a time-domain signal and A = fft(a), np.abs(A) is its amplitude spectrum and np.abs(A)**2 is its power spectrum. The phase spectrum is obtained by np.angle(A).

Then I took a look at DCT, but I could not figure out how to derive the phase from the results, because only real values are contained, so the np.angle() will be $0$. From the FFT results I can take out interesting elements, calculate amplitude, phase and plot them as cosines:

[A*np.cos(2*np.pi*f*t + ph) for t in range(N)]

where A is the amplitude, f is the frequency, ph is the phase-shift and t is the variable, i.e. time. So my question is: can the same be derived from DCT results and how?

Not really, as the transform is real. However, one could interpret the sign as a poor man's phase, being "quantized" or restricted to values $0$ or $\pi$. In other words, $1 = 1.e^{0.\imath}$ and $-1 = 1.e^{\pi.\imath}$.

[EDIT] There are some instances where people use the sign of DCT or (real) wavelet coefficients, for subpixel image registration or matching, as alluded to in DCT Sign Only Correlation and Its Application to Image Registration:

  In the DCT case, all DCT coefficients are real values, as a result,
  the phase information corresponds to the sign of DCT coefficients.


Colored noises: Stationary or non-stationary? I know for sure that white noise is considered as a stationary sound, but is that true for the rest of them? http://en.wikipedia.org/wiki/Colors_of_noise
Also, can we meet some of these sounds in nature? Is there any relation between these noises and noises that you can hear in a town(vehices,traffic,people chatting) or are these noises just an artificial invention?

Thanks in advance for your ideas and thoughts!:)

Stationary processes have a spectrum that's time-independent. As the Wikipedia link shows, the colored noises each have a characteristic spectrum, are therefore not time-dependent, and thus stationary.


Hamming Window LUT function I'm working on a project where I need to generate a LUT for Hamming Window to perform a FFT.

I'm using this function (doing the table with MATLAB):

n = 0:4095;
LUT_length = length(n);
WIN_HAM_v = round(4095*(0.54-(0.46*cos(2*pi*n/(LUT_length-1)))));

I have some doubts if the function is correct

Your code is correct.

This can also be implemented by calling the built-in MATLAB function hamming:

WIN_HAM_v = round(4095*hamming(LUT_length));

You can (should) compare the two methods yourself:

win1 = 4095*(0.54-(0.46*cos(2*pi*n/(LUT_length-1))));
win2 = 4095*hamming(LUT_length);
maxError = max(abs(win1(:)-win2(:)));
% result: 2.2737e-12


How does gaussian convution stacking work? If I first convolve an image with a Gaussian filter with σ = 1, and then convolve
the output of that image with a Gaussian with σ = 2

Does this give the same effect as convolving the original image with σ = 3?

From this Wikipedia page, you can get:

So for your example, $\sigma_X = 1$ and $\sigma_Y = 2$. The result has a standard deviation of $\sigma_Z = \sqrt{\sigma^2_X + \sigma^2_Y} = \sqrt{1 + 4} = \sqrt{5}$.


Hermitian symmetry in OFDM systems I am trying to understand the usage of Hermitian symmetry in OFDM systems and have a couple of questions regarding this.

What is the reason of using the Hermitian symmetry in OFDM?
How can we arrange the data in terms of Hermitian symmetry?
In OFDM system where do we arrange our data as Hermitian symmetry?
Is there any reference or documentation that explains Hermitian symmetry in OFDM system?

The Hermitian symmetry is used to obtain a real-valued time-domain signal. It is a special case of OFDM called discrete multitone (DMT). It exploits a property of the discrete Fourier transform (DFT), namely that the DFT of a real-valued signal has Hermitian symmetry. The motivation is usually the channel: if the signal shall be transmitted over a low-pass channel without additional radio-frequency upconversion (which requires an oscillator at the transmitter and a local oscillator at the receiver) the transmit signal must be real-valued.

Given an OFDM system with $N$ subcarriers the following conditions must hold
$$
c_{-k} = c^*_k\\
c_0\in \mathbb{R}
$$
where $c_k$ is the complex value of subcarrier with index $k\in[-N/2\ldots N/2-1]$. This can be achieved by mapping $N/2$ complex symbols to subcarriers $0$ to $N/2-1$ and then assigning the respective complex conjugate value to the according subcarrier with negative index. This would normally take place after the mapper and before the inverse DFT (IDFT) block.


FFT of sine signal with phase discontinuities I am detecting sine signals within a specific bandwidth [0,B] Hz,

The received sine have multiple phase discontinuities, as shown below

I want to measure the frequency of the sine signal but I get a distorted spectrum as expected

What's the best practice to get rid of the high leakage? I have tried windowing but did not really work, I also tried FIR low-pass filter with f_c = B Hz, but also did not work.

Since the time series appears to have a very high SNR and you have a good idea of the frequency, know the sample rate, and can detect the discontinuities,

$$ \min_{A,f,\theta} \left\lVert y_i - A \sin(2\pi\frac{f}{f_s} i + \theta)\right\rVert^2_2 \quad \{i\in \text{outside discontiuties}\}$$

will be locally convex if the minimization is near the true values.  The DFT is good at pulling tones out of noise. You don't seem to have a lot of noise so, the DFT isn't always the best tool.

not enough room in comments

If you really want to use a DFT based technique on segments, try Macleod's frequency interpolator:

M. D. Macleod, "Fast nearly ML estimation of the parameters of real or complex single tones or resolved multiple tones," in IEEE Transactions on Signal Processing, vol. 46, no. 1, pp. 141-148, Jan 1998.
doi: 10.1109/78.651200
Abstract: This paper presents new computationally efficient algorithms for estimating the parameters (frequency, amplitude, and phase) of one or more real tones (sinusoids) or complex tones (cisoids) in noise from a block of N uniformly spaced samples. The first algorithm is an interpolator that uses the peak sample in the discrete Fourier spectrum (DFS) of the data and its two neighbors. We derive Cramer-Rao bounds (CRBs) for such interpolators and show that they are very close to the CRB's for the maximum likelihood (ML) estimator. The new algorithm almost reaches these bounds. A second algorithm uses the five DFS samples centered on the peak to produce estimates even closer to ML. Enhancements are presented that maintain nearly ML performance for small values of N. For multiple complex tones with frequency separations of at least 4π/N rad/sample, unbiased estimates are obtained by incorporating the new single-tone estimators into an iterative “cyclic descent” algorithm, which is a computationally cheap nonlinear optimization. Single or multiple real tones are handled in the same way. The new algorithms are immune to nonzero mean signals and (provided N is large) remain near-optimal in colored and non-Gaussian noise
keywords: {Fourier analysis;Gaussian noise;amplitude estimation;frequency estimation;harmonic analysis;interpolation;iterative methods;maximum likelihood estimation;optimisation;phase estimation;signal sampling;spectral analysis;white noise;AWGN;Cramer-Rao bounds;FFT;amplitude;cisoids;colored noise;complex single tones;discrete Fourier spectrum;fast nearly ML estimation;frequency;frequency separations;interpolator;iterative cyclic descent algorithm;noise;nonGaussian noise;nonlinear optimization;nonzero mean signals;parameter estimation;peak sample;phase;real single tones;resolved multiple tones;sinusoids;unbiased estimates;uniformly spaced samples;Amplitude estimation;Computational efficiency;Frequency estimation;Iterative algorithms;Maximum likelihood estimation;Noise level;Parameter estimation;Phase estimation;Phase noise;Signal processing algorithms},
URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=651200&isnumber=14197


PyWavelets CWT: normalization? Vs Scipy? Related. The equation being implemented normalizes by sqrt(1 / scale):
$$
C_{a, b} = \frac{1}{\sqrt{a}} \sum_k s(k)\left( \int_{-\infty}^{k+1} \overline{\psi \left(\frac{t-b}{a} \right) }dt - \int_{-\infty}^{k} \overline{\psi \left(\frac{t-b}{a} \right) }dt \right)
$$
so why does pywt._cwt do * sqrt(scale)? Also, how is it normalizing (L1, L2, etc)?

Edit: After reviewing Scipy's cwt, I realized it's implementing exactly what I described as "ideal" in the answer, thus it's same as comparing the two (except a few obvious details). I don't know in which scenarios pywt's algorithm is better, but generally - easy win for Scipy. I've opened a PR to address pywt's key flaw.
Edit2: Found an apparent edge of pywt's over scipy's in the least expected place: low scales. Addressed here.
Edit 3: after fixing the low precision in pywt, I actually found pywt superior to scipy in all of a variety of synthetic signals, but squeezepy's always beat both. Tests here, with more coming.

Let's follow the math from incubation to delivery. It begins with psi, a rescaled morlet2 (as shown previously) at a scale $a=64$, and $\sigma=5$:
$$
\psi = \psi_{\sigma}(t/a) = e^{j\sigma (t/a)} e^{-(t/a)^2/2} \tag{2}
$$
gets integrated and L1-normalized: -- (see caveat2 below)
$$
\psi = \psi_{\text{int}}(t) = \frac{1}{a} \int \psi_{\sigma}(t/a)\ dt \tag{3}
$$
The norm is off by a scalar factor, but what matters is it's same for different $a$; you can verify it's so (for large $a$, increase integration limits). Note that it's the wavelet that's L1-normalized, not the integrated wavelet $\psi_\text{int}(t)$; the distinction is essential, as we'll see.
$\psi$ is then  convolved with the signal $s$, and finite-differenced and negated; from here, this is:
$$
 -\text{DIFF}(\text{conv}) = [s \cdot (\psi_{\text{int}}(t - 0) - \psi_{\text{int}}(t - 1)),\ ...] \tag{4}
$$
Inspecting closely, this is equivalent to undoing the former integration, with certain numeric differences (addressed later), so we're left with
$$
 \psi = \frac{1}{a} \psi_{\sigma}(t/a) \tag{5}
$$
Lastly, * sqrt(scale):
$$
 \boxed{\psi = \frac{1}{\sqrt{a}} \psi_{\sigma}(t/a)} \tag{6}
$$
which is thus ultimately L2-normalized. We can verify this on a test signal; L2 is biased in amplitude/power representation toward lower frequencies, suppressing high frequencies. If above is right, not only can we see this, but apply a fix by removing * sqrt(scale) from code:

Implementation-specific caveats discussed next.

How reliable is the normalization across scales? -- Not great (-- code2):

On right-hand side we see the best case using limited number of samples, on left we see the resampling implementation. Clearly due to inherently poor resolution of small number of samples, our norms don't do well for low scales - but the recomputed variant at least maintains a mean norm approximately constant across them.
More interestingly, the resampled one suddenly appears quadratic after $\log_2(a)=6$; what's the deal? Recall "It begins ... at scale $a=64$"; this is when we resample with more samples than in original (int_psi):

whose diff is

In 1LR we see int_psi at scale=256 is resampling the same values from int_psi multiple times, which manifests as zeros in diff (2LR).
Next we square the values to get an idea of the L2 norm, and compare against the L2-normalized values in 3LR; we find normed values blowing unnormed values out of water. Is this excessive, or appropriate?
Compare against recomputed at same scale; in 4LR, we at least see recomputed's shape - it's no longer negligible. This suggests the resampling to be overestimating wavelet amplitude, and since in norm computation this is significantly offset by half the values being zero due to diff, then the overestimation must be great.

How great? We can compare peak wavelet amplitudes across scales of the unnormed resampled diff(int_psi_256) and recomputed w (morlet2/scale); norming changes both by the same factor (*sqrt(scale)), won't make a difference.

Surely enough, again at $a=64$ we see critical behavior; peak amplitude of resampled saturates at a minimum, whereas the recomputed continues to decrease. The difference at $a=256$ is more than great enough (~12db=4x, if squared -> 16x) to offset the zeros in diff, hence why the L2 norm grows.
But notice the L1 norm remains well-behaved at higher scales. To find out why, we must precisely compare the proportion of zeros with the overestimation factor. We already know latter to be 4x; for former, I took a non-flat region (otherwise it's zeros anyway) and found ratio of zeros to number of samples: 0.75. Which means, 1 in 4 values are non-zero, and each are 4x greater ... the exact proportion needed to retain the L1 norm.
That zeros proportion isn't very clear from above plots - so instead, dots:

If diff undoes cumsum, why do either at all?
Good question. In fact, by ridding of both, I found the quality degradation (zipper) effect disappear:

But something else also disappears - the high frequency scales (left), and the sign flips (I undid this by removing minus at coef = -). I partly offset this by replacing * sqrt(scale) by 1 / sqrt(scale) (right), but that gets us to the L2 normed case (1 / scale is overkill), and is a baseless modification in light of everything found so far. So can it be done? Sure, one can implement without the int-diff dance, but I won't be exploring it here.

Caveat: effective norm
Recall, conv gets trimmed to len(x), so the entire wavelet may not participate in computation, whereas above computations assume otherwise. The assumption is truer the longer x is - but here's what happens if len(x) == 512, with same wavelet precision=10 (2^10=1024 int_psi), compared against recomputed wavelet:

A mess! Both distort significantly, possibly implying it's inherent to great enough scales pinned against short signals. A simple workaround is zero-padding the signal, but I wouldn't vouch for it without understanding the caveats (which I don't for wavelets, but do for DFT).
Furthermore, note that the norm isn't only per-scale, but also per-shift ($b$); less of the wavelet is in the signal when it's just entering as opposed to when its center is at signal's center.
But keeping it per-scale, I confess I cheated in the earlier example with $f=4$; initially I picked $f=1$, but the L1 norm plot didn't look nice, and I didn't understand why. Since now I do, here are both L1 and L2 per scale for $f=1,64$:

Compare against earlier plot; note how the left modes (low scales) of both norms are unaffected, only the high scale (low frequency) ones are now much greater (L2) and lower (L1), just as expected.

Caveat 2: * step doesn't do the L1-normalizing
Instead, it results from diffing the resampled integrated wavelet; the * step only sets the correct initial scale, and can be omitted without losing the norm. This does not happen with recomputing the wavelet, or resampling the unintegrated wavelet.
(Note that recomputed scipy.signal.morlet2 is L2-normalized, but only because it's explicitly normalized; the diff-resample-int is "self-normalized").


how to determine a length of a channel? ( given path, delay) Given a channel matrices for each subcarrier of a frequency selective channel (OFDM): i.i.d..
channel =(randn(ntx,nrx, N)+1i*randn(ntx,nrx, N))/sqrt(2); simulated via channel simulator and a number of paths, delay

With these parameters, Can I define a length of the channel or a number of clusters $L$?

With these parameters, Can I define a length of the channel

There's now $n_{tx} \times n_{rx}$ (or, depending on how you consider the this, $n_{rx}$), with potentially different lengths, not just one channel!

The only thing that makes sense to define as "the channel length" is the maximum length among these many channels.

However:

The fact that you're using a working OFDM already says someone else has found the length of maximum length of the channel, and used that as length of the cyclic prefix.

So, really, nothing to do here. If there was anything to do, your OFDM would either be inefficient (too much CP) or wouldn't work (too short CP).


what is smoothing in very basic terms What is smoothing and how can I do it?

I have an array in Matlab which is the magnitude spectrum of a speech signal (the magnitude of 128 points of FFT). How do I smooth this using a moving average? From what I understand, I should take a window size of a certain number of elements, take average, and this becomes the new 1st element. Then shift the window to the right by one element, take average which becomes the 2nd element, and so on. Is that really how it works? I am not sure myself since if I do that, in my final result I will have less than 128 elements. So how does it work and how does it help to smooth the data points? Or is there any other way I can do smoothing of data?

EDIT: link to followup question

Smoothing can be done in many ways, but in very basic and general terms it means that you even out a signal, by mixing its elements with their neighbors. You smear/blur the signal a bit in order to get rid of noise. For example, a very simple smoothing technique would be, to recalculate every signal element f(t) to as 0.8 of the original value, plus 0.1 of each of its neighbors:

f'(t) = 0.1*f(t-1) + 0.8*f(t) + 0.1*f(t+1)

Note how the multiplication factors, or weights, add up to one. So if the signal is fairly constant, smoothing doesn't change it much. But if the signal contained a sudden jerky change, then the contribution from its neighbors will help to clear up that noise a bit.

The weights you use in this recalculation function can be called a kernel. A one-dimensional Gaussian function or any other basic kernel should do in your case.

Nice example of one particular kind of smoothing:

Above: unsmoothed signal
Below: smoothed signal

Examples of a few kernels:


Determining time variance of integrator Say we have an integrator $y(t)= \int_{- \infty}^{t} x(t) dt$

It is incorrect to use the following method to determine if it is time variant/invariant?

First, when we have input signal $x(t),$

$y_1(t)= \int_{- \infty}^{t} x(t) dt$

Then, we shift $y(t)$ by $t_0$

$y_1(t-t_0)=\int_{- \infty}^{t-t_0} x(\tau-t_0) d\tau$

If we have input signal $x(t-t_0),$ we get

$y_2(t)=\int_{- \infty}^{t} x(\tau-t_0) d\tau$

And since $y_1(t-t_0)\neq y_2(t)$

Then the system must be time variant.

This is not correct. First, you have to write it like that:

$$
\begin{align}
y(t) = \int_{-\infty}^t x(\tau)d\tau
\end{align}
$$

I.e. your variable in the integrand is different from the integral limits. Accordingly, your calculation will be different, as you will only replace $t$ with $t-t_0$ and not $\tau$.


Is there a spectrum sensing library? I was building my cognitive engine for a Cognitive Radio and, first things first, I have to make my CR aware of the environment by sensing the spectrum.

I read about Matched Filter Detection, Energy Detection and Cyclostationary Feature Detection: while I found plenty of libraries for machine learning, I did not found any implementation of the three sensing techniques mentioned above.

So, is there a spectrum sensing library (C++ would be nice)?

There's a whole load of spectral estimation libraries available for GNU Radio. I'm personally fond of gr-specest, the GNU Radio spectral estimation toolbox.

It does a lot popular algorithms, inluding Cyclostationary spectrum estimators, ESPRIT, MUSIC, Welch, etc.

For "energy detection" and "matched filtering", you wouldn't need any add-on library – GNU Radio already comes with the necessary blocks (here: complex to magnitude square, and FIR filters) in highly optimized versions.

GNU Radio is written in C++, but you can connect blocks using Python, if you prefer so, and write your own blocks in either python or C++.

I'd recommend going through the official GNU Radio Tutorials, as they will pretty quickly give you a feeling how you can implement the sensing for a CE with GR.

In fact, if you search  for (Cognitive Radio) AND (GNU Radio) on IEEExplore, you'll find quite some papers (634 at the time of writing). It might be a good idea to investigate those since 2014 to get a feeling for the state of the art.


Identifying Bad Sectors on magnetic stripe, or noise in signal I am starting with an audio signal picked up from the microphone jack of my iPhone, generated by running a standard magnetic strip ATM card through a sensor.

On an old card, there will be a few degenerate bits:

Now filtering out this damage is one issue, which I have asked here: Removing noise from F2F signal

But how to determine the integrity of a particular section of track?   How could I mark a particular area as being damaged?

I basically need some measure of smoothness for any given sample.  My gut feeling is to work with derivatives.  Taking the derivative amplifies noise.

But I can't see it clearly. Can anyone?

Typically you wouldn't directly use time-domain samples of a signal to declare its integrity (or lack thereof). I'm not familiar with the magnetic-stripe specification, but in most systems, any encoded data is protected by some form of checksum. This checksum is usually located in the same stream of symbols as the data of interest. Many systems also include error-correction coding, to provide tolerance to some number of symbol errors while still being able to recover the encoded information successfully.

Once the signal has been demodulated and you have recovered the symbol stream, you would calculate the checksum of the recovered data and compare it to the checksum that you read from the medium itself. If they match, then you can assume that the data is correct. If they do not, then you can try re-reading the medium again, or after some number of tries, you might declare it as a bad or damaged region.

This all begs the question, though: what is your goal of marking these regions as damaged? What will you subsequently do with the damaged areas?


Segment highly packed nuclei with dotted high intensity signal inside each object I collected high magnification images on cell nuclei and I need to segment them (red circle in the image). The edges of the nuclei are visible but inside each nuclear object there are bright spots of positive signal (Not an imaging artefact, if I will retake the images it will be the same).

I tested different way to segment the image (using python scripts based on scikit-image, scipy.ndimage and opencv):
- Threshold->morphology (close/opening)-> watershed (after gaussian filtering or median filtering to try to remove the dots)
- Graph-cuts (Al-Kofahi algorithm implemented in the Farsight package: https://github.com/gbkedar/farsight)

All the methods lead to oversegmentation because of the internal bright spots.

To solve that I thought was a good idea to isolate the dots with a laplacian of gaussian and used the opencv function cv.inpaint to fill the dots with the intensity of the surrounded pixel values. The result was not good either.

The questions are:
- Do you think that the best strategy is to isolate the bright dots and then try to replace them with the surrounding signal?
- Any strategy that will lead to a better result?

Any help will be greatly appreciated

I wouldn't bother with correct bright pixel segmentation.
I think such a strategy would add complexity and another stage in the algorithm which might add errors at the eventual product.

I would use some form of scaling function for the pixel value so that the bright pixels would blend-in better while maintaining contrast in the darker pixels. One such function could be the log function:

  im = log(double(im));

Another one might be the square root.
You might even get away with using a max value for the image , e.g.

  thr = 150; im(im> thr) = thr;

or a combination of both, for instance a capped transform function
thr = 100;
tform_im = img.^0.5;
tform_im_scaled = tform_im - min(tform_im(img>thr)) + thr;
img(img>thr) = tform_im_scaled(img>thr);

Either way good luck!


Weighted Sum of Auto Correlation - Lower Bound Given a vector $ v $ with elements $ {\left\{ {v}_{n} \right\}}_{n = - \infty}^{\infty} $ and denoting $ {v}_{n}^{\left( k \right)} = {v}_{n - k} $, namely, a shifted vector by $ k $ elements (Mind the vector is infinitely long).

How could one prove that there exist $ \alpha > 0 $ such that:

$$ \sum_{k = -\infty}^{\infty} {2}^{- \left| k \right| } \left \langle {v}^{\left( 0 \right)}, {v}^{\left( k \right)} \right \rangle \geq \alpha {\left\| v \right\|}^{2} $$

One could see it a weighted sum of the Auto Correlation Function of the vector $ v $.

Hence can be written as:

$$ \sum_{k = -\infty}^{\infty} {2}^{- \left| k \right| } \sum_{n = -\infty}^{\infty} {v}_{n} {v}_{n - k} \geq \alpha {\left\| v \right\|}^{2} $$

My understanding of the question is that there some bounderay on how slow the Auto Correlation decays (There is nothing special about $ {2}^{-\left| k \right|} $ there, it can be replaced by any $ \frac{1}{c} > 1 $).

Thank You.

Defining $ a \left[ k \right] = {2}^{- \left| k \right|} $.
Moreover, the Auto Correlation function of $ v $ defined as $ {r}_{vv} \left[ k \right] = \left \langle {v}^{\left( 0 \right)}, {v}^{\left( k \right)} \right \rangle = \sum_{n = -\infty}^{\infty} {v}_{n} {v}_{n - k} $.
Pay attention that Auto Correlation is Hermitian Function.

Using the definition of Convolution one could write:

$$ \left( {r}_{vv} \ast a \right) \left[ 0 \right] = \sum_{k = -\infty}^{\infty} {2}^{- \left| k \right| } \left \langle {v}^{\left( 0 \right)}, {v}^{\left( k \right)} \right \rangle $$

Using the Convolution Theorem one could write that:

$$ \left( {r}_{vv} \ast a \right) \left[ 0 \right] = \int_{- \pi}^{\pi} {R}_{vv} \left( \omega \right) A \left( \omega \right) d \omega $$

Where $ R \left( \omega \right) $ and $ {R}_{vv} \left( \omega \right) $ are the DTFT of $ {r}_{vv} \left[ k \right] $ and $ a \left[ k \right] $ respectively.

One should notice the DTFT of $ a \left[ k \right] $ is defined only one sided. Yet since its symmetrical it can well calculated:

$$
\begin{align*}
A \left( \omega \right) & = DTFT \left\{ a \left[ k \right]  \right\} = \sum_{k = -\infty}^{\infty} a \left[ k \right] {e}^{-j \omega k} = \sum_{k = 0}^{\infty} {2}^{-k} {e}^{-j \omega k} + \sum_{k = 0}^{\infty} {2}^{-k} {e}^{j \omega k} - 1 \\
& = \frac{1}{1 - 0.5 {e}^{-j \omega}} + \frac{1}{1 - 0.5 {e}^{j \omega}} - 1 = \frac{1 - {c}^{2}}{1 - 2 c \cos \left( \omega \right) + {c}^{2}} = \alpha > 0 \quad \forall c < 1
\end{align*}
$$

In the above $ c = {2}^{-1} = 0.5 $ yet actually this will hold for any $ c < 1 $.

So the integral is given by:

$$
\begin{align*}
\int_{- \pi}^{\pi} {R}_{vv} \left( \omega \right) A \left( \omega \right) d \omega & = \int_{- \pi}^{\pi} {R}_{vv} \left( \omega \right) \frac{1 - {c}^{2}}{1 - 2 \alpha \cos \left( \omega \right) + {c}^{2}} d \omega \\
& \geq \alpha \int_{- \pi}^{\pi} {R}_{vv} \left( \omega \right) d \omega = \alpha {\left\| v \right\|}^{2}
\end{align*}
$$

As requested.

By the way the result must be real since $ a \left[ k \right] $ is symmetric and $ {r}_{vv} $ is hermitian function and hence its transform is real.


Inverse system of a system with integral I am trying to find the inverse system of the following (I tried finding the mathematical inverse function but since it is not the same I am not so sure) . Can someone help me find it?

$$ y(t)=\int_{t-1}^{t+1}\cos\left(\frac{\pi\tau}{8}\right)x(\tau)d\tau $$

It seems that we have the same homework. You probably are Greek.

After having the same question, I came to the conclusion that in order for this system to be invertible, you have to prove that for any given $$x_1(t),  \, x_2(t)$$ the following sentence is correct: $$ x_1(t) \ne x_2(t) \,\, {\Rightarrow} \,\, y_1(t) \ne y_2(t) $$

As like you proved when a math function is $"1-1"$.

For this very example, if you have the following the following two signals:

$$
\
  x_1(t) = δ(t-4)
\
$$

and

$$ x_2(t)=0 $$

it turns out that while initially $ x_1(t) \ne x_2(t) $, we have $ cos(πt/8)x_1(t) = cos(πt/8)x_2(t) $, for any $t$. This means that

$$ \int_{t-1}^{t+1}\cos\left(\frac{\pi\tau}{8}\right)x_1(\tau)d\tau = \int_{t-1}^{t+1}\cos\left(\frac{\pi\tau}{8}\right)x_2(\tau)d\tau $$

and the system is not invertible.


How can a RLS algorithm utilise Wiener filter as FIR (M-tap) block? I'm currently working with a dataset of $5000$ pulses of $N=15000$ samples each. I managed to implement the RLS algorithms with a FIR M-Tap filter such that $M\leq 15000$ ($150$ seems to achieve the best results).

I would like to initialise the coefficients of my filter to the values achieved via Wiener (optimal filter of $N$ coefficients).

Please see link here on how to implement the Wiener Filter and a simple explanation of its operation in terms of deconvolution (described as an LMS filter using the Wiener-Hopf equations, but is a Wiener Filter implementation. It should be noted that the true Wiener filter is an MMSE solution which the LMS solution will converge to).

Compensating Loudspeaker frequency response in an audio signal


Noise removal using adaptive noise cancellation algorithms in real time systems I am newbie in signal processing. I want to remove noise from real time audio signal using adaptive noise cancellation(LMS&RLS) algorithms. I don't know how to choose reference signal(noise signal) for the algorithm from the incoming noise corrupted signal.  So can anyone help me to select the noise signal for the reference of ANC algorithms?
Thanks in advance

In real time application adaptive noise cancellation uses single sensor for signal reception.The single sensor consist of a  delay  z-Δ to produce a delayed version of incoming signal d(n), denoted by x(n), which de-correlates the noise while leaving the target signal component correlated. So here not needed reference signal separately but delay creates uncorrelated noise from the incoming noise corrupted signal.


Confusion calculating length of Kaiser window I am reading Digital Signal Processing Using Matlab by Ingle and Proakis (3rd ed.)
In Chapter 7, Section 7.3,
I am confused why he is adding 1 to main lobe width for calculating length of window.
I have attached snapshots below.
First/above snap is of matlab code of example 7.8 and 2nd snap is of code of example, both snaps contain the confusing line as highlighted.

This is just an empirical formula found by Kaiser for determining the necessary filter length for a given transition width. That formula is given as Equation $(7.30)$ on page $332$:
$$M=\frac{A_s-7.95}{2.285\,\Delta\omega}+1\tag{1}$$
I think that Kaiser came up with a formula for determining the filter order (hence without the $+1$ in the equation), and the authors of your book preferred to have a formula for the filter length, so they took the original formula and added $1$ to it.
Judging from some of your previous questions, you seem to be confused when it comes to the terms filter order and filter length. For FIR filters, filter length is the number of coefficients (taps). Filter order is the (minimum) number of delay elements necessary to implement the filter. It's just like with polynomials: their order is one less than their number of coefficients. E.g., a second-order polynomial has $3$ coefficients:
$$P_2(x)=a_2x^2+a_1x+a_0\tag{2}$$
Coming back to FIR filters, you always have
$$\textrm{filter length}=\textrm{filter order}+1\tag{3}$$


Best Car Recognition Algorithm Could anyone suggest the best algorithm for a real-time Car recognition (say in a parking space)?
I am planning to implement the same on FPGA as well. Kindly suggest.
Thank you.


How to calculate pitch fundamental frequency f( 0) ) in time domain? I am new in DSP, trying to calculate fundamental frequency ( f(0) ) for each segmented frame of the audio file. The methods of F0 estimation can be divided into three categories:

based on temporal dynamics of the signal time-domain;
based on the frequency structure frequency-domain, and
hybrid methods.

Most of the examples are estimating fundamental frequency based on the frequency structure frequency-domain, I am looking for based on temporal dynamics of the signal time-domain.

This article provides some information but I am still not clear how to calculate it in the time domain?

https://gist.github.com/endolith/255291

This is the code, I have found, used so far :

def freq_from_autocorr(sig, fs):
    """
    Estimate frequency using autocorrelation
    """
    # Calculate autocorrelation and throw away the negative lags
    corr = correlate(sig, sig, mode='full')
    corr = corr[len(corr)//2:]

    # Find the first low point
    d = diff(corr)
    start = nonzero(d > 0)[0][0]

    # Find the next peak after the low point (other than 0 lag).  This bit is
    # not reliable for long signals, due to the desired peak occurring between
    # samples, and other peaks appearing higher.
    # Should use a weighting function to de-emphasize the peaks at longer lags.
    peak = argmax(corr[start:]) + start
    px, py = parabolic(corr, peak)

    return fs / px

How to estimate in time domain?

Thanks in advance!

This paper by Philip McLeod details finding pitch in the time domain, and specifically that it avoids the reduced resolution of equivalent frequency domain approaches due to the windowing required to eliminate spectral leakage.

https://cdn.instructables.com/ORIG/F3J/3Q7U/J20P2R5T/F3J3Q7UJ20P2R5T.pdf

I don't have experience with this approach but would be very interested in comparative results to the best frequency domain approach- if you could post them, assuming you pursue this further.


applying single pole filter to simple signal How to apply single pole filter to simple signal

fs = 500;
ts = 1/fs;
t = 0:ts:1-ts;
fc = 50;
S = cos(2*pi*fc*t);

If your pole is at $z=p$ (real or complex), this will do:

Y = filter( [1], [1, -p], S );

Free bonus: if you also want a zero at $z = c$:

Y = filter( [1, -c], [1, -p], S);

This last expression corresponds to the filter:

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


what is the relation between a minimum phase system and a FIR system? I want to see whats the relation between being an FIR system and a minimum phase system and when an FIR system can be minimum phase too with respect to have linear phase

The only FIR system that is both minimum phase and linear phase is where the output is the same as the input: y[i] = x[i]


Advantage of complex filtering in multirate applications I've seen it mentioned in passing in various papers on signal processing and filter design that complex FOR filters can be more efficient when it comes to multirate applications. However I cannot think of a scenario where this is applies. Of the systems I've seen, the baseband is complex which means generally the baseband for interpolation (for example) filters makes most sense with real coefficients so +ve and -ve frequencies are both passed.

Can anyone give examples where a complex bandpass filter makes sense, and specifically to areas like decimation and interpolation.

For one practical example, I'll point to GNURadio's frequency translating FIR filter block: https://github.com/gnuradio/gnuradio/blob/master/gr-filter/lib/freq_xlating_fir_filter_XXX_impl.cc.t

It's a channelizer  & decimator block that gains an efficiency by spinning a user provided FIR LPF up to be a complex BPF where the channel is and filtering first vs. spinning the channel down first.  The efficiency comes from performing the frequency rotation at the lower sample rate after filtering and decimation vs. before.


how could window function really suppress those frequencies that we don't want? I know we can use window function to decrease leakage phenomenon by making the end of a non-periodic signal to be zero at both end, so that to make it like "continuous and periodic" when doing fft.

However, when we multiply window function with original signal
(ex: w(t)xs(t)), it's like doing convolution (w(f)*s(f)) in frequency domain, just like w(f) function slide from left to right and integrate with the area with s(f).

And my question is that since it would still integrate the area which is not the real frequencies of the original signal, and have values after convolution, how could this method really suppress those frequencies?

What make me confused is because I think using window function is the same as the method of doing FIR, which also multiple a window function with the ideal low pass signal in frequency domain. And the FIR characteristics are related to how many areas during integration (ex: the size of transition band).

Thank you very much!

You are right that the effect of windowing is the same for the design of FIR filters and for spectral analysis. You will always have a compromise between main lobe width and relative side lobe level (relative to the main lobe). The main lobe width determines the frequency resolution, and the relative side lobe level determines the amount of spectral leakage. A rectangular window has a narrow main lobe (desirable), and relatively large side lobes (undesirable). Other windows achieve different compromises between those two parameters. In FIR filter design, the equivalents to frequency resolution (main lobe width) and leakage (side lobe level) are transition band width and stop band rejection.

You should read the classic paper by Fred Harris: On the Use of Windows for Harmonic Analysis with the Discrete Fourier Transform.


Fourier Transform of a Pulse Train (Frequency Comb Simulation) I am trying to take a Fourier transform of a femtosecond pulse train and plot the frequency spectrum.
I followed this question and am able to create a pulse train and plot both the time and frequency domain.
My question is, why do the pulses overlap in the time domain graph, since their width is on the order of femtoseconds, and the distance between them (T) is much larger?
Is it an issue with my time vector or am I interpreting the graph incorrectly?

%Using a Gaussian carrier envelope for the e-field, I show a comb of laser
%frequencies defining frep, and its representation in freq space. The
%sampling frequency is 10 times the freq (f=10*fs, fs=fr so fr=f/10).

%Define variables
fr=120e6;             %rep rate, units: Hz
fs=10*fr;             %sampling frequency, units: Hz
T=1/fs;               %period, units: s (~fs) (sampling time)
tau=100e-15;           %proportional to pulse duration, units: s
%tau=(pulse_duration)/sqrt(2*log(2));
wl=2*pi*fr;           %carrier frequency, units: Hz
ceo=0;                %offset freq., can change later

%Define electric field
L=1500;               %sets axis length
t=(0:L-1)*T;          %time, units: s

Es=zeros(length(t));
for n=0:100           %n is the number of pulses (depends on L)
    for i=1:length(t)
        Es(n+1,i)=exp(-(t(i)-n*T)^2/tau^2)*exp(-1i*wl*t(i));
    end
end
E=sum(Es,1);

%Fourier Transform to frequency
F=abs(fft(E));        %Fourier transform the vector Es
FT=F(L/2:L-1);           %Take out half due to symmetry
FTn=FT/max(FT);        %Normalized y-axis
f=(1/L)*(L/2:L-1)*fs;    %Frequency for plotting

%Plots versus time and frequency
figure(1)
plot(t,real(E))
xlim([0 max(t)/2])
xlabel('Time (s)')
ylim([-1.2 1.2])
ylabel('Amplitude')
figure(2)
plot(f,FTn,'k')
xlabel('Frequency (Hz)')
xlim([min(f) max(f)])
ylabel('Normalized Amplitude')

Let us first see the theory.  By pulse train I guess you mean a periodic function from minus infinity to plus infinity.  Your pulse train will have the Fourier series(FS) in stead of Fourier transform(FT).  Fortunately the FT of a single pulse and the FS of the pulse train are related.  Let the pulse train be periodic with P Hz.  Then the FS of the pulse train will be sampled version of FT with a sampling rate of P Hz.  i.e. except for a scaling factor, the FS values are the samples of the FT sampled at P Hz. And the single pulse FT is a sinc function.  Thus pulse train will have FS coefficients which are samples of the sinc function.


How do I calculate the range of a fixed-point number with $a$ integer bits and $b $ fractional bits? What does fixed-point number range represents? Why we use formula $2^a - 2^{-b}$, why minus $2^{-b}$? Where $a$ is number of integer bits, and $b$ is number of fractional bits.

If we have for example $a = 8$ and $b = 2$, don't we have a possibility to represent $2^a + 2^{-b}$ number and so that will represent the range?

Let's assume that we are dealing with unsigned number types. If you would use all $a+b$ bits for the integer part then the set of possible numbers would be:
$$\left\{0, 1, 2, 3, \dots, 2^{a+b}-1\right\}.$$

These numbers can be divided by $2^b$ (or multiplied by $2^{-b}$) to take use of $b$ bits for the fractional part, resulting in this set of possible numbers:

$$\left\{0, 1\times2^{-b}, 2\times2^{-b}, 3\times2^{-b},\dots, \left(2^{a+b}-1\right)\times2^{-b}\right\}
\\=\left\{0, 2^{-b}, 2\times2^{-b}, 3\times2^{-b},\dots, \underline{\underline{2^a-2^{-b}}}\right\}.$$

So your formula gives the largest number that can be represented (double underlined).

The formula can also be understood as going a step $2^{-b}$ or one least significant bit (LSB) worth backwards from $2^a$ which is the first number that has the same truncated binary string representation as 0, the first number in the system. In a similar way in the 8-bit unsigned integer system we take a LSB-sized step backwards from 256 (1 0000 0000 binary) to obtain 255 (1111 1111 binary) which is the largest representable number in that system.

I just found out that range has an established meaning in arithmetic: Largest value minus smallest value. The smallest value happens to be zero so in this case the largest value equals the range.


Find tracing/plotting stacked bands/lines in an image I have images with lines/bands, the bands/lines are parallel but can be of different gray levels, I want to extract or trace the lines/bands as single line, all lines/bands are parallel.

Source:

Output Plot:

The source image contain sine wave lines (10 pixels in amplitude) that are stacked on top on of each other forming what looks like bands. I used the location of the max of the cross-correlation of the first column with other columns of the image.

I am looking for other approaches as cross correlation is too expensive and very sensitive to noise.

Why not try just summing the columns of the image?

Well, you could try smoothing the image first, but that will have the same effect as just smoothing the output.

R code for this is below the first lot.

Here's the extra image, processed.

R Code Below

# 26593
# install.packages("png")

file <- "Q26593_original.png"

library(png)

img <- readPNG(file)

d <- dim(img)

par(mfrow=c(2,1))
plot(1:2, type='n', main="Plotting Over an Image", xlab="x", ylab="y")
lim <- par()
rasterImage(img, lim$usr[1], lim$usr[3], lim$usr[2], lim$usr[4])

plot(colSums(img[,,1]), type="l")

R Code #2 Below

# 26593
# install.packages("png")
# install.packages("EBImage") --- Didn't work.
# ## try http if https is not available
# source("http://bioconductor.org/biocLite.R")
# biocLite("EBImage")

file <- "Q26593_original.png"

library(png)
library(EBImage)

img <- readPNG(file)

d <- dim(img)

#xx <- (0:(d[2]-1) - d[2]/2)/d[2]
#yy <- (0:(d[1]-1) - d[1]/2)/d[1]
#sigma <- 0.2

#w <-  exp(-yy*yy/sigma) %*%  t(exp(-xx*xx/sigma))
#w1 <- exp(-yy*yy/sigma) %*% t(rep(1,d[2]))

f <- makeBrush(11, shape='disc', step=FALSE)
imgf <- filter2(img, f)
imgf <- imgf/max(imgf)

par(mfrow=c(3,1))
plot(1:2, type='n', main="Filtered Image", xlab="x", ylab="y")
lim <- par()
rasterImage(imgf, lim$usr[1], lim$usr[3], lim$usr[2], lim$usr[4])

#plot(1:2, type='n', main="Window Function", xlab="x", ylab="y")
#lim <- par()
#rasterImage(w1, lim$usr[1], lim$usr[3], lim$usr[2], lim$usr[4])

plot(colSums(imgf[,,1]), type="l", main="Column sum of filtered image")

plot(filter(colSums(img[,,1]), rep(1,11)), type="l", main="Column sum of unfiltered image, filtered")


what is the C/D equation in dsp? I have one example test question by this kind:

So in this ebook (Digital Signal Processing Using MATLAB) as the DSP references and not found C/D block equations on it?

and this c/d search results?

Thanks for your attention.

C/D is simply a "continuous-to-discrete" time conversion. In other words, sampling! :-) You can use $T$ to sample your continuous time signal $x(t)$ into a discrete time signal $x[n]$. Just replace $t=nT$ in your $x(t)$ equation, to turn it into a discrete time signal $x[n]$. D/C is the opposite operation. I think you can move on from here. Let us know if you can work it out now or if you need further help.


Understanding FFT: FFT size and bins I'm trying to implement an FFT to understand how it works.
I'm using a random sound file (so i can't control the size of the signal), with a sampling frequency of 44,1Khz.
I'm using this library

  https://www.nayuki.io/res/free-small-fft-in-multiple-languages/Fft.java

The objective is to apply this formula to get the frequency:
F = n * Fs/N

With n number of bins.
Fs sampling frequency, and N size of the FFT

1)An Fft is supposed to have a length, most of them use a power of 2 radix. But how can i know the length of the FFT if i apply it to an entire array of data ? is it the size of that array ?

2)I don't understand the bins, how is it processed ? i know that FFT works with chunks of data, but if i give a full array of data, how are the bins "created".

Question 1

If you apply it over the entire length of the array, the length of the FFT would be the length of the array. But, the FFT is more efficient if the length is a power of two, so it is common to pad 0's onto the end of the signal until its length is a power of 2.

Overly simple example...

x = [3.4, 2.56, 1.3]

x has a length of 3, the next power of 2 after 3 is 4, so we change x to be

x = [3.4, 2.56, 1.3, 0]

and apply an FFT with length 4.

Another big BUT! If your signal is long, it becomes extremely inefficient to do the whole thing at once. You would not want to try to do an FFT on an audio file the length of even a short song. In that case, we break the signal into chunks of some reasonable size, perform an FFT on each, and average the results.

Odds are good that what you actually want to do with your data is a not just a standard FFT, but rather the averaging process I described above. Google Bartlett and Welch methods for more details.

Question 2

I'm not 100% sure what you're asking about here. I'm going to interpret it as you wanting know how the width of the frequency bins are determined and run with that.

The width of each frequency bin is determines solely by the rate the signal was sampled at and the length of the FFT. The width of each bin is the sampling frequency divided by the number of samples in your FFT.

df = fs / N

Frequency bins start from -fs/2 and go up to fs/2.
That means if sampled at 100Hz for 100 samples, your frequency bins will be width 1Hz. If you take 200 samples, you will now have 2x as many frequency bins and their width will be 1/2Hz each.