DSP Concepts Visually Explained Similar to this question: Visually stunning math concepts which are easy to explain, what are some great visualizations of basic DSP concepts such as FFTs, filters, etc?

I don't know if it qualifies as quite "visually stunning", but you might want to check out my blog article:  DFT Graphical Interpretation: Centroids of Weighted Roots of Unity

The concept of the $1/N$ normalized DFT as a center of mass calculation was a major "aha moment" for me.  It is a good answer for "What does the DFT really mean?"

By request, here is one of the figures from my article:

A little explanation is in order.  The top graph is a time domain representation and the polar graphs on the bottom are on the complex plane.  The left most circle is for bin zero, aka the DC bin, the second is bin one, and so on.  The little blue circle is the center of mass and is also the bin value as as complex number.

$$ A = \frac{2 \pi n}{N} $$

This sample has 3 cycles per frame with a phase of 3.  Bin three (the fourth polar graph) clearly shows the bin value has magnitude of $1/2$ and the phase value of 3 is almost $\pi$ and therefore almost halfway around the circle.

There are many more examples, and more thorough explanations with the math in the article.


gabor wavelet and gabor transform I was reading a paper 'A Flexible Framework for Local Phase Coherence Computation' (paper URL) on using wavelets for local phase coherence. I just want to be sure. Is the gabor wavelet the same as the gabor function in the gabor transform (STFT)? If not, can we use a Windowed Fourier Transform instead?

The Gabor wavelet is (almost) the same as the Morlet wavelet (some authors distinguish these two by an additional constant to fulfill the wavelet's admissibility condition). And yes, it is the same function as in the Gabor transform (gaussian windowed oscillation).

The application of phase coherence should be independent of the underlying wavelet and the nature of the transform (Wavelet, Gabor/STFT, other weird metrics). Essentially important is the concept of reproducing kernels and the intrinsic (Heisenberg) uncertainty that influence you measure of phase coherence and must be taken care of.


What does it mean that a zero is slow? I am studying control systems, and have encountered the definition of a slow zero. I am searching on internet and in books, but I don't understand the meaning this. I know that if a zero is too slow, it introduces a resonance peak and an overshoot in the step response, but what is a slow zero?
I thought it was a zero close to the imaginary axis in the root locus, so a zero at low frequencies, but I am not sure, since I cannot find a definition.
I know that a pole that is close from the imaginary axis settles quicker than a pole far from the imaginary axis. But I did some simulations in Matlab, and using a lead compensator, I have seen that if I decrease frequency of the zero in the lead compensator, the overshoot in the step response increases.

Here, the lead compensators I have used are as below.
lead = (1+10*0.05*s)/((1+0.05*s));          %red line
lead_2 = (1+30*0.083*s)/((1+0.083*s));      %blu line
lead_3 = (1+50*0.083*s)/((1+0.083*s));      %green line

Can somebody help me?

Again I refer you to the Murray's book. Which is free by the way

http://www.cds.caltech.edu/~murray/books/AM08/pdf/fbs-limits_18Aug2019.pdf

From the book "Furthermore a zero is said to be “slow” if its magnitude is smaller
than the intended closed loop bandwidth."


Karhunen loeve transform question I have read some about Karhunen-Loeve Transform (KLT) and its application to the field of seismic data processing.
The method as I understand it based on decomposing the data (actually mostly used in image processing) using SVD on its covariance matrix and projecting the data back after manipulating one or some of the covariance matrix's eigenvalues.
One of the assumptions, when using on seismic data to remove random noise, is that the features to keep, should be aligned in time. If this is not the case, then one have to correct for the observed move-out.

For ex. in the top picture below 400 time samples (y-axis) for more than 50 data traces(x-axis) are shown. Each trace represents data recorded on one geophone. There are multiple "events" in the data represented by horizontal lines across the figure. There is no time delay for the observed event across the traces.
The picture in the bottom is the result after applying KLT on the data, by zeroing the smallest eigenvalue of data's covariance matrix.
My question is, why the assumption of time alignment of the "event" (horizontal features)? The method fails if the events occurred along a line with a slope, that is the same event was observed with a delay from one trace to the next.

Thanks

There is something that is not clear of what you have done with the data, and that is who do you form the random vectors to perform de SVD (or EVD) on the covariance matrix.

1 -The KLT can be succesfully used on a one dimensional signal (only one Geophone), taking frames of $M$ samples and estimating a covariance matrix from it, and the performing Eigenvalue decomposition, or with the SVD of the Data matrix (which supposedly is the same), then you sort the eigenvalues (or singular values) in descending order, and remove the smallest ones according to a particular criteria. With this you will remove part of the noise.
Based on the process i described, there is no assumption about any correlation between the signaals on the geophones (and in reality i think there will be, but i am not an expert in geology), so if you perform the same processing to every signal, it doesnt matter what delay you have between particular events, the process will still remove noise.

2 - But, if you construct your random vector with one sample from each geophone taken at the same time instant, then the time delay might affect the correlation between each sample from a geophone and then applying KLT might not reduce the noise as it is expected.

The approach number 1 i described is called Temporal PCD (Principal component analysis) and the second one is called Spatial PCA, assuming you put the geophones in different locations.

Personally i think it should not fail if you do spatial PCA, because even though there will be time delay, the samples from different geophones at same time instants will have some degree of correlation, except if you put them hundreds of miles apart, but if they are all in the same region, then they will each sense approximately the same signal but contaminated with noise and probably convoluted with an impulse response of the geographical region, but as long as there is correlation, you can apply KLT.


Yes an FFT approach will work to identify frequency content, but you may need longer data sets to see something useful. For the FFT you will see the frequencies from 0 to your sampling rate (25000), and there will be the same number of samples in frequency as you have in time.  Therefore your frequency resolution is 25000/N where N is the number of samples.

Note, importantly the frequency resolution is dependent only on the time length of your data, T and is independent of the sampling rate. The above relationship works out as if you sample more, N will be proportionally larger in a fixed block of time T.

Therefore

$$ Res = 1/T $$

Where Res is the frequency resolution of your measurement and T is the time length of the data.

For your two datasets, you have 122 data samples for the first one, and and 222 samples for the second; which is 4.42 ms and 8.84 ms of data. The frequency resolution is 25000/122 = 204.92 Hz for the first block data and 25000/222 = 112.61 Hz for the second block (or 1/4.42 ms and 1/8.84 ms).

Let's look closer at your data, the second set that was larger. From a plot in the time domain (see below), we see a very large transient at the initial part of your capture that would likely obscure our ability to see the general frequency content, so I chose to work with a subset of your data from 2.8ms to 8.84 ms (which would decrease the frequency resolution further unfortunately but allows us to not be overcome by the initial impulse).

Showing the region after the initial impulse:

Here is the FFT of that portion of your sequence. I had removed the mean before taking the FFT which is why there appears to be a large null in the center.

The peak is near 500 Hz, with around 170 Hz of resolution.

This is what the FFT looks like without removing DC or the initial transient in the beginning.

Windowing your data first can help reduce the side-lobes in frequency to give you more dynamic range (if you are looking for lower level signals that may be masked by the side-lobes from larger signals), but know that that will decrease your frequency resolution further.

Here is the same plot as above of your entire second data set after applying a Kaiser window.

And a zoom in of above showing the FFT bin sample locations:

Zooming in further, we see the detail below. The Kaiser window used was with a Beta of 8, which has a frequency resolution that is approximately 3x wider than the non-windowed example. Since I used the complete data set (and the effects of the impulse are masking the lowest frequency region accordingly), the resolution non windowed would be 112.6 Hz, and with the window used around 330 Hz.

My recommendation to get more meaningful results is to increase the time length of your capture (assuming your frequency is stationary for a longer time). Given the frequency roll-off visible from these FFT's, it is not necessary to sample at 25 KHz; so you could reduce your sampling rate if that helps you increase the time length of the capture.


Is it possible to implement FIR band pass filter using two FIR band stop filters? I am trying to filter a signal using a given FIR band stop filter. I have to implement a band pass filter.  Assuming I want to get the signal between normalized cutoff frequencies [0.2,0.3] - can I use two band stop filters over the same signal - one for the range [0,0.2) and second for the range (0.3,1] ?

Yes, there's no reason why you couldn't do that. If you cascade two systems in series, the overall transfer function (and therefore frequency response) is the product of the two systems' individual responses. Therefore, you can apply the first filter to attenuate frequencies in the range $[0,0.2\pi]$, then the second will attenuate frequencies in the range $[0.3\pi, \pi]$. The result is that all frequencies except for those in the range $(0.2\pi, 0.3\pi)$ are attenuated, which looks a lot like a bandpass filter response.

Note: Don't take the choices of brackets/parentheses in the intervals above too seriously. Any realizable filter has a finite-width transition region, not instant transitions between passband and stopband.


using Short-Time Energy and zero crossing rate I am trying to identify a sound region of a word from a wav file (44100Hz). My research gave me the idea that I must use "Short-Time Energy" and "zero crossing rate". But the problem is I am not clear with the formulas to code. So if someone can guide me to a much clear Pseudo code I will greatly appreciate it. And also if someone can tell me the order of using this I will highly appreciate that too. Thanks in advance.

Audio processing FFT of chunks of audio
This will give you the short term energy as a function of time and frequency.

The zero crossing in itself is easy to code (the example is conceptual, but very inefficient)

x0 = x(1);
for i = 2:length(x)
  x1 = x(i);
  if (x0 > 0 && x1 < 0) || (x0 < 0 && x1 > 0)
     % you have a zero crossing
     doSomethingWithAZeroCrossing(i);
  end
  x0 = x1;
end

However, zero crossings are fairly sensitive to benign and typically irrelevant changes in the audio signal, so zero crossings have only limited usefulness. In particular

even small amounts of noise can generate lots of extra zero
crossings and "jitter" the existing ones
Fluctuations in DC offset or very low frequency noise (like an Air conditioner or a heating fan picked up my a microphone) move the zero crossings all over the place


Audio sounds oddly distorted https://clips.twitch.tv/TangentialObliqueLobsterSaltBae
https://clips.twitch.tv/CharmingBashfulAlfalfaRickroll
https://clips.twitch.tv/NaiveKindGorillaGivePLZ

(Disclaimer: These were not recorded on or uploaded from my computer.)

In the 1st video linked above, the water sloshing sounds heavily distorted and the same type of "noise" heard throughout that segment of the clip can be heard periodically throughout this person's videos, including the other ones I linked.

I've been trying to help this person hear the problem with his audio, but he's not hearing it and I'm not entirely sure what the cause of it could be or how one could even describe it.

EDIT: He switched to different broadcasting software and finally hears how bad it sounds now. Now we just need to know the cause.

It's mono, not stereo like one would expect, and it sounds like exacerbated compression artifacts to me, like in the difference signal between the stereo channels of a lossy compressed audio, but I'm not sure if that explains everything.

Given that he uses an external mixer for recording, there may be a simple explanation: He has connected a stereo cable to a balanced input of the mixer. The balanced input will typically take the left channel as "hot" and right channel as "cold", and will effectively see the mono difference signal: hot minus cold. He should use a stereo to 2x mono splitter cable, or use a stereo input in the mixer if available.


Understanding Communication Theory Jargon for Statistical Linear Models Coming from a statistical background, I understand the following formula:

\begin{equation}
y = X\beta + \epsilon
\end{equation}

as a linear regression model, with data input-output pair: $(X,y)$, noise term, $\epsilon$, and regression coefficient vector $\beta$.

Now I am confused in how this system is described in communication theory. I see terms such as Rayleigh Fading, MIMO, Multi-access, Multi-antenna, Band-limited etc ... in relation to equations that look exactly like statistical linear regression.

So far I suppose I can work out that for communication theory,

\begin{equation}
y = Hx + \epsilon
\end{equation}

$x$ is some input signal at a transmitter antenna, $\epsilon$ is noise again, $y$ is the signal at the receiver antenna ... but what is $H$? I don't understand how this term should be described, and how it is modelled with respect to all those previous terms I have listed? Reading journal papers I see the simple system as above described with lots of jargon I can't follow, and I'm unaware of the implications this has on the modeling.

If someone could clarify all this jargon and connections for me that would be suuuper helpful!

The equation $y =  hs + n$ is a discrete-time model of a certain kind of wireless communications system.

$s$ is a complex number, drawn from a finite set called a "constellation", that represents the information being transmitted.
$n$ is a complex, circularly-symmetric Gaussian random variable with zero mean and variance usually denoted by $N_0$.
$h$ represents the effect of the communications channel. In this simple model of a wireless channel, $h$ is a complex Gaussian random variable with zero mean and variance one. This model is called Rayleigh-fading because $|h|^2$ is a Rayleigh-distributed random variable.
$y$ is the received complex number. From it, the receiver must make an estimate of $s$. This estimate is usually calculated by first estimating $h$, and then finding $$ y' = \frac{h^* y}{|h|^2} = s + \frac{h^*n}{|h|^2} $$ Here you can see the importance of $|h|^2$ to enhance the noise, which is why the channel is called Rayleigh. The estimate of $s$ is calculated by finding the constellation element that is closest (in Euclidean-distance sense) to $y'$.


How much to decimate? What I don't understand is the connection between the sample rate, the baud rate, decimation and using convolution or correlation.

I'm sampling at 2048000 Hz but I don't know my baud rate, I estimated it (by looking at the time length of one of the short pulses) to be about 2000 pulses per second.

So these are my steps (roughly):

get samples
demodulate
decimate by about 1000
correlate with ones

Here is an example of what my signal looks like after demodulation:

You can see that there are 4 symbols but many many samples.  How can I turn those many many sample into just 4 symbols? I think it is related to convolution/correlation, but I'm not sure how.

Here's a general example on how to perform the convolution + decimation for an ASK signal.

amplitudes = [1, 0 0 1 0 1 1 0 1 0];
N = 100; % the symbol length.
pulse = 1.2*ones(N, 1); % the Rx pulse (used for correlation)

% the tx pulse (used to recreate your RX signal);
tx_pulse = pulse + 0.2*sin(2*pi*10/N*(1:N)');

% artificially create a signal
signal = zeros(N*length(amplitudes),1);
for i=1:length(amplitudes)
    signal((i-1)*N+(1:N)) = amplitudes(i) * tx_pulse;
end
signal = signal + 0.01*randn(size(signal));

subplot(3,1,1);
plot(signal);
title('RX signal');

subplot(3,1,2);
after_conv = conv(signal, pulse);
plot(after_conv);
title('After convolution');

subplot(3,1,3);
stem(after_conv(1:N:end));
title('After decimation');


How to find the value of regularization parameter/s when having one or multiple priors I am working on a image reconstruction framework, specifically interferometry. In simple words we have data that is pertubated by noise, say

$$
V^o_i = V_i + \sigma_i
$$

Where $V^o$ is a complex value and usually called the measured/observed visibility which is corrupted by Gaussian noise, $V$ is the Fourier Transform of an image $I$ and $\sigma$ is the noise of each measurement.

We can solve the problem solving a non-linear equation such as:

$$
 \Phi(I)= \arg \min_{I} ||\frac{V^o - V^m(I)}{\sigma}||^2
$$

where we call $V^m$ model visibilities. I don't want to extend much more here, but since we have an infinity number of solutions and given that a chi-squared will return a noisy image we want an image that represents the highest possible fidelity to the data.

In that sense we want to add prior to regularize the solution. Here, we can think of Maximum Entropy, L1 Norm, Total Squared Variation, etc. Consequently, the equation turns into:

$$
 \Phi(I, \lambda)= \arg \min_{I} ||\frac{V^o - V^m(I)}{\sigma}||^2 + \lambda P(I)
$$

where $\lambda$ is the regularization parameter and P is the prior function.

I've been reading papers about how to find the $\lambda$ parameter and most of them use an estimation calculating the variance of the noise. But in that case then $\lambda = \text{Var}(\sigma)$? This is also called the discrepancy principle (See Regularization in Image Restoration and Reconstruction by W. Clem Karl) In that case, $\lambda$ would be the same for each prior right? But it has been seen doing the L-curve plots for same dataset and different prior that the "best" $\lambda$s for each one of the cases are different.

Moreover, since I am constructing a framework, the idea is to plug different priors, and then the equation will turn to:

$$
 \Phi(I, \lambda)= \arg \min_{I} ||\frac{V^o - V^m(I)}{\sigma}||^2 + \lambda_1 P_1(I) + \lambda_2 P_2(I)+..+\lambda_n P_n(I)
$$

On that case the calculation of the different $\lambda$s is not that trivial as before, I think....

I would appreciate if you can give me a slight hand with both of these cases.

Cheers

I have just found this paper that calculates the $\lambda_i$ parameters on an L-surface: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.135.5969&rep=rep1&type=pdf


Partial derivative of image using mask I want to find the partial derivative of an image with respect to its X coordinate. I read that it can be done by convolution of the image with a mask like:

             fx = conv2(im1, 0.25 * [-1 1; -1 1])

I do not understand how is it calculating the partial derivative - can anyone explain?

Partial derivatives in discrete domains are computed using differences (finite difference in more proper terms). There are many methods to do that such as central finite difference or forward, backward and so on.

The convolution operation you mentioned computes forward distances in horizontal direction and results in x-derivatives (namely the $x$ component of the gradient).

To also compute in $y$ direction you might do:

fx = conv2(im1,0.25* [-1 -1; 1 1])

0.25 is there just to scale the values to appropriate range (weights each pixel the same).

For more information on finite differencing you might check:
http://en.wikipedia.org/wiki/Finite_difference#Forward.2C_backward.2C_and_central_differences

Convolution is a local-masking operation where at each pixel you compute the dot product of the kernel and matrix part that corresponds to the kernel. Then, kernel is shifted. I won't go to the details of convolution as it is already explaind in http://en.wikipedia.org/wiki/Convolution. However, let's say we have a matrix

$
   I=
  \left[ {\begin{array}{cc}
   1 & 2 & 4 \\
   6 & 14 & 17\\
   16 & 41 & 18\\
  \end{array} } \right]
$

If we convolve this with the kernel mentioned and discard for now the boundary conditions, we end up with

$
   G=
  \left[ {\begin{array}{cc}
   -2.25 &  -1.25 \\
   -8.25 &  5.0 \\
  \end{array} } \right]
$

(this is part in the middle of the result).

So, you should notice that basically, per each element, we compute a sum of 2 neighbors vertically on the left, sum of 2 neighbors vertically on the right and subtract these intermediate results. As subtraction is done horizontally, it computes the horizontal derivatives. This is what convolution intrinsically does.


Shift a signal by a given number of samples Is shifting a signal by a given number of samples, say $n$ samples, the same as the following in matlab: signal(n:end) instead of signal(1:end).

That is, if the size of $signal$ is 200x1, does it make sense to remove the first $n$ rows and consider the size of $signal$ at (200-n)x1 for the shifted signal?

No, this is not shifting but windowing by a rectangular window. Shifting will not change the values of the samples, it will just change the relation between the sample values and the corresponding time indices.

In Matlab a (finite length, discrete-time) signal is simply represented by a vector of numbers. There is no relation to time indices. The interpretation which sample value corresponds to the index $n=0$ is up to you. So shifting a signal without any further processing in Matlab (or any other software) is meaningless.

If it's just about plotting, you can use stem(n,x) where x is your signal vector, and n is the corresponding vector of time indices. Shifting just involves specifying a different vector n, the data in x remain of course unchanged.

Shifting becomes important if it is part of some processing involving more than one signal, e.g. the computation of the deterministic auto-correlation or cross-correlation. In such a case you simply have to implement the corresponding formula.

EDIT:

If you have two signal vectors x and y which have the same number of elements N, and you want to add a shifted (e.g., delayed) version of y to the signal x then you can do two things. If y is shifted by M samples, you could define the sum by

z = [x(:),zeros(M,1)]+[zeros(M,1),y];

In this case the indices where only one of the two signals in non-zero (due to the shifting operation) are also taken into account. This may or may not be useful. If you want to disregard these regions and only take into account indices where both signals overlap, you can do this:

z = x[M+1:N]+y[1:N-M];

This is simply the middle region of the other sum vector z defined above.


Output of a band pass filter Signal

x(t) = | cos(2πft)|, f != 0

is the input of an ideal band-pass filter  with bandwidth B = 3f/2.
How can i Find the Output of the y(t) of the signal.

Your signal is periodic with period $\frac{1}{2f}$ hence has a Fourier series expansion with spacing between the frequency components of $2f$. The output of your ideal filter then is the sum of the Fourier components that fall within the bandwidth.

Because your bandwidth is less than the spacing between the frequency components your output will contain at most one frequency component.


the rate of noise presented in images I suppose that an image is affected by a multiplicative noise (or additive noise and etc.). I want to quantify and measure the rate of noise (for exemple 1%, 2%, ...,25%) present in an image. I used the snr and psnr but they give me the values measured in dB (decibel). Also, how can I implement this in MATLAB?

Can I use the mse or the average between the noisy image and the original image mean?

After the answer of Mr. Tolga Birdal, I put the following code in matlab:

I=rgb2gray(imread('lena.png'));
J=imnoise(I,'speckle', 0.1); % 0.1 is the variance of the noise
Sigma1 = estimate_noise(I)

Its give to me :

  Sigma1 = 6.9360

Sigma2= estimate_noise(J)

and

  Sigma2= 10.1409

this mean that the level of noise present in the image J equals to 10.14 % ?

You could check out my MATLAB function, which estimates the variance of the noise. It is also mentioned in a Mathworks blog post, in comparison with other methods.

Regarding $\sigma$, you mention in the question, simple answer is yes, longer answer is no. It refers to the variance of the noise normalized with image dimensions as follows:

$$\sigma = \frac{\sigma_0\sqrt{0.5\pi}}{6(W-2)(H-2)}$$

$-2$ is due to the convolution kernel, applied previously. However, this is not a final quantity, giving you the absolute noise level. Nevertheless, for the variance, you could use it as a percentage, just like you mentioned.


Causality and ROC of a stable LTI system So I am looking at a stable LTI system whose input is $x[n]$ and output is $y[n]$. The equation relating the two is here:
$$
y[n-1]-\frac{10}{3}y[n]+y[n+1]=x[n]
$$

I was able to compute its system function $H(z)$ to be:
$$
H(z)=\frac{z^{-1}}{z^{-2}-\frac{10}{3}z^{-1}+1}
$$

I then performed partial fraction expansion to end up with:
$$
H(z)=\frac{\frac 38}{\left(1-3z^{-1}\right)} - \frac{\frac 38}{\left(1-\frac 13z^{-1}\right)}
$$

What I am struggling with is I can't decide what its ROC should be since the system is specified to be stable which means it has to include the $j\omega$-axis. I also need to determine if the system is causal or not. Since a system can only be causal and stable when all its poles are on the left hand plane, I am not sure.

Any help determining both the ROC and causality of the system would be very helpful.

This looks like homework, so I'll give you a hint to help you figure this out by yourself.

There are 2 poles: one at $z=3$ (i.e., outside the unit circle) and one at $z=\frac13$ (i.e., inside the unit circle). For such a system there are three possible ROCs; one corresponds to an anti-causal system, another one to a causal system, and a third one to a non-causal system with a two-sided impulse response. But only one of those ROCs includes the unit circle, and that's the ROC corresponding to a stable system.


Get frequency and amplitude from an audio sample I'm wondering how can i get the frequency and the amplitude from a 16 bits sample.
I have actually an array of 16 bits signed int and I would like to get the frequency and the amplitude of each sample. I've heard from FTT and stuff but I'm kinda new to signal processing about audio.
Any advice would be appreciated,

Thanks.

I'm wondering how can i get the frequency and the amplitude from a 16 bits sample

That makes no sense. A single sample is just a number.

A number doesn't have a frequency.

The number itself is the amplitude.


Variance Due to white noise input I have the problem below. It sounds simple but for some reason I have been stuck on it for a long time and don't know what am doing wrong.

Am trying to solve this using correlations.
So we all know that the variance $\sigma _y^2\:=\:R_y\left(0\right)$ which is basically the correlation of $y$ evaluated at zero.
Now we know the convolution $y(n)\:=\:x(n)* h(n) = \sum _{m\:=-\infty }^{\infty }h\left(m\right)x\left(n-m\right)$
Using the idea of correlation, I get the following:
$R_y\left(k\right)\:=\:E\left(y\left(k\right)y\left(n+k\right)\right)\:=\:E\left[\sum _{m=-\infty }^{\infty }h\left(m\right)x\left(k-m\right)\sum _{m=-\infty }^{\infty }h\left(m\right)x\left(n+k-m\right)\right]$
Now am so stuck here and just don't know what to do. Please help me by directing me in the right direction and by providing clear and coherent steps and explanations. Thank you very much in advance.

You have a typo in your definition of $R_y(k)$, and an error in the time indices when developing the equation. The proper definition of autocorrelation for complex signals (in the case of wide-sense stationary processes) is
$$R_y(k)=\mathbb{E}\left[y(n)\overline{y(n+k)}\right],$$
and setting $k=0$, we obtain
$$\begin{align}
R_y(0)&=\mathbb{E}\left[y(n)\overline{y(n+0)}\right]=\mathbb{E}\left[\left|y(n)\right|^2\right]\\
&=\mathbb{E}\left[\left|\sum_{m=-\infty}^{\infty}h(m)x(n-m)\right|^2\right]\\
&=\mathbb{E}\left[\left(\sum_{m=-\infty}^{\infty}h(m)x(n-m)\right)\overline{\left(\sum_{i=-\infty}^{\infty}h(i)x(n-i)\right)}\right]\\
&=\mathbb{E}\left[\sum_{m=-\infty}^{\infty}\sum_{i=-\infty}^{\infty}h(m)\overline{h(i)}x(n-m)\overline{x(n-i)}\right]\\
&=\sum_{m=-\infty}^{\infty}\sum_{i=-\infty}^{\infty}h(m)\overline{h(i)}\,\mathbb{E}\left[x(n-m)\overline{x(n-i)}\right]\\
&=\sum_{m=-\infty}^{\infty}\sum_{i=-\infty}^{\infty}h(m)\overline{h(i)}\,R_x(m-i)\\
\end{align}$$
where I first use the fact that the product of sums is the sum of products (and I use a different summation index to better show the cross-products), and second, that the impulse response of the system is deterministic and therefore constant for the purposes of taking expectation.
Now, we can group those terms where $m=i$ and those where $m\neq i$, obtaining
$$\begin{align}
R_y(0)&=\sum_{m=-\infty}^{\infty}\sum_{i=-\infty}^{\infty}h(m)\overline{h(i)}\,R_x(m-i)\\
&=\sum_{m=-\infty}^{\infty}h(m)\overline{h(m)}\,R_x(0)+\sum_{m=-\infty}^{\infty}\sum_{i=-\infty\\ i\neq m}^{\infty}h(m)\overline{h(i)}\,R_x(m-i)\\
&=R_x(0)\sum_{m=-\infty}^{\infty}|h(m)|^2+\sum_{m=-\infty}^{\infty}\sum_{i=-\infty\\ i\neq m}^{\infty}h(m)\overline{h(i)}\,R_x(m-i)\\
\end{align}$$
Finally, we know that the input noise is white, and therefore $R_x(k)=0$ for $k\neq0$. We also know that $R_x(0)=\sigma_x^2$. Thus, the second term vanishes, and we obtain
$$\begin{align}
R_y(0)
&=R_x(0)\sum_{m=-\infty}^{\infty}|h(m)|^2=\sigma_x^2\sum_{m=-\infty}^{\infty}|h(m)|^2
\end{align}.$$


Design of a digital A-weighting filter with arbitrary sample rate I want to A-weight a time series with arbitrary sample rate.

An analog A-weighting filter is defined exactly by IEC 61672-1. But there's no definition for a digital filter. One method is to use the bilinear transform (BLT) to convert the analog filter to the digital filter (as done here Applying A-weighting). However this method suffers from extreme warping near nyquist (even when the analog poles/zeros are pre-warped):

Figure 1: A-weighting frequency response comparison where the sample rate is $25600\textrm{ Hz}$.

Instead I'm thinking of using an algorithm than can design a digital IIR filter with arbitrary frequency response and plugging in the frequency response of the analog A-weighting filter.

Is this a good approach?
If so, is there a particular algorithm that would be well suited for this?

I've looked into MATLAB's yulewalk but I would need a corresponding Python implementation to try out. I've also come across Berchin's FDLS method in a few places, like this question for instance, but all of the links appear to be broken.

It's a common misconception that the approximation of an analog filter by a digital filter must be bad close to Nyquist. This idea might come from the ubiquity of the bilinear transform, for which this is usually indeed the case. Of course, there are certain constraints on the frequency response of discrete-time filters at Nyquist, but they do not necessarily need to result in a bad approximation of an analog filter in that frequency range. The quality of the approximation close to Nyquist depends on several factors, among which are the properties of the frequency response of the analog filter, and the fact whether only the magnitude or also the phase of the analog filter need to be approximated.

I've designed a 6th order IIR filter approximating the frequency response of an analog A-weighting filter, as defined here:

$$H(s)=\frac{k\cdot s^4}{(s+129.4)^2(s+676.7)(s+4636)(s+76655)^2}\tag{1}$$

with $k=7.39705×10^9$.

I chose a sampling frequency of $48$ kHz. The design procedure is a heuristic iterative procedure I came up with some time ago. It's a least squares approximation based on the equation error method, and I might write up all the details some day.

Below is a plot of the design result. Note that both plots go up to Nyquist ($24$ kHz). The left figure shows that one can't see any difference between the logarithmic plots of the analog and the digital frequency responses. The right-hand figure shows the approximation error, defined as the absolute value of the difference of the magnitude responses. The error shows a typical least squares behavior.

You can check out the filter yourself. Here are the coefficients:

b =

   0.169994948147430
   0.280415310498794
  -1.120574766348363
   0.131562559965936
   0.974153561246036
  -0.282740857326553
  -0.152810756202003

a =

   1.00000000000000000
  -2.12979364760736134
   0.42996125885751674
   1.62132698199721426
  -0.96669962900852902
   0.00121015844426781
   0.04400300696788968


Find half-angle of complex exponential Suppose we have a complex signal $\cos(2\theta) + j\sin(2\theta) = e^{j2\theta}$. Is there any way to find half-angle $\cos(\theta) + j\sin(\theta) = e^{j\theta}$ complex signal without going into polar form and preferably with no divisions? Iterative algorithms are fine.

For small $\Theta$ an approximation would be to simply halve the imaginary value.

You could also use the ratio of the real and imaginary values as the index to a lookup table that would return the half angle value.  You could not store every possible value, of course, so you would have to choose some set of values and then lookup the value that is closest to the actual ratio.


Sox, mixing mono inputs into mono output (efficiently) I want to mix two mono files, signal.wav and noise.wav, into a single mono file, signal_noise_mono.wav.

Currently, I know how to mix the two input files into a stereo file and then remix it to mono, but I would like to mix down to mono directly, and avoid the intermediate file.

# Merge signal (mono, L) and noise (mono, R) into a stereo file
sox -M signal.wav noise.wav signal_noise_stereo.wav

# Remix channels 1 & 2 from the stereo file into a single mono file
sox signal_noise_stereo.wav signal_noise_mono.wav remix 1,2

Any tips? Have I missed an obvious sox switch?

As far as I can understand, you want to add two signals, together?
I guess that you are missing -m switch.

Generate two test files (sweep and sine):

sox -Dr 16000 -c 1 -b 16 -n sine.wav synth 5 sin 1000
sox -Dr 16000 -c 1 -b 16 -n sweep.wav synth 10 sine 100:8000

Merge them together:

sox -m sweep.wav sine.wav -Dr 16000 out.wav

Using soxi out.wav you can see that produced file is mono:

Input File     : 'out.wav'
Channels       : 1
Sample Rate    : 16000
Precision      : 16-bit
Duration       : 00:00:10.00 = 160000 samples ~ 750 CDDA sectors
File Size      : 320k
Bit Rate       : 256k
Sample Encoding: 16-bit Signed Integer PCM

Now you can observe the spectrograms:

sox sine.wav -n spectrogram -o sine.png
sox sweep.wav -n spectrogram -o sweep.png
sox out.wav -n spectrogram -o out.png

Sine:

Sweep:

Merged together:


What is meant by "information encoded in time domain"? Author on pg 292 claims the windowed-sinc "is not a filter for signals with information encoded in time domain", because its "step response has overshoot and ringing". That phrase, along "encoded in frequency domain", are used in text, but I'm unsure they're ever explicitly defined. My guess:

Frequency-domain encoded: periodic processes in a system - something that recurs, persists the duration of the signal, and is sensibly described by Fourier coefficients
Time-domain encoded: the opposite - time-local events in a system, such as a sharp rise, temporary change in DC offset, momentary noise injection.

A single square pulse would qualify for latter, as the physical system that generated it almost certainly didn't synch hundreds of sinusoidal sources just right. A noiseless EEG would qualfiy for former, as neuronal excitations recorded as polarity alignments have periodic traits.
Is this accurate? If so, how is a filter's step response relevant in describing time-local events - and are there other "responses" characterizing time-encoded information?

Note: I'm aware of the step response's role in control systems, circuit design; I ask specifically about non-periodic "events" in a signal (if that's what time-encoded means).

The simplest (imho) explanation is this.
Consider a time-domain signal which has narrow pulses or sudden jumps in value or on/off switches at certain (unknown!) time instants. Such a signal is said to carry its information in those specific time instants at which such sudden changes occur.
A PWM (pulse width modulation) signal is an excellent example for such class of signals.
If you consider the frequency domain reflection of those sudden changes in the signal amplitudes, you will see that they have wide band spectrums. They are local in time, therefore spread in frequency.
When you filter such a signal with a typical LPF (low pas filter) or alike, which discards parts of the spectrum away and passes only a subset of it, you may effectively be throwing away part of the information that's carried by the signal too.
For example, the specific instants at which the PWM wave switches on/off becomes blurry, and hence the resolution of the information is reduced.
The author seems to mention this issue with filtering of time-domain encoded information...


Designing the filtering for a time domain data acquisition application I want to filter sensor signals for a data acquisition application on an ATmega microcontroller. The interesting parts of the signal are in the range 0-5Hz, and I'd like to leave these mostly unchanged, but strongly attenuate 50Hz power line hum and anything above that.

The signal will later be evaluated in the time domain, especially the time and amplitude of large peaks (typically lasting well over a second), and the area under those peaks.

Since the signal is generated by a different device with its own power supply, I expect to see some power supply hum as well as general noise.

In order to minimize hardware components I plan to use a simple RC filter (-3dB at 18Hz) for antialiasing, combined with a fairly high sample rate (> 1kHz), and then apply digital filtering. The signal will then be sent on to a PC at a much lower sampling rate (probably <10Hz). However, I'm not very experienced with filtering, so I'm a bit overwhelmed with the choices available.

Given the fast sampling rate and very limited CPU power, filters that require many floating point operations would at best be a good deal of work to implement. I played around with filter design tools a bit and looked at the properties of Butterworth and Bessel filters or just chaining several first-order filters in a row (since they can be implemented very efficiently), and finally found an interesting section about the humble "moving average" filter on dspguide.com.

After playing around in Octave a bit, it seems to me that by chaining four moving average filters, two with 1/50s and two with 1/60s window, would give me

around 0.5dB attenuation at 5Hz, the outer limit of what might still be interesting signal
over 100dB attenuation at both 50Hz and 60Hz, even if the frequency is 1% off
over 80dB attenuation for anything above 100Hz
linear phase shift in the signal band, for a fixed delay of <40ms

This looks very promising, but since I only started playing with all of this very recently, I'm not completely sure my results are valid. Is this a good filter for my task, or would you recommend something else?

There might also still be noise left in the range of 5Hz to 50Hz which would cause aliasing when the signal is sampled at 10Hz by the PC. Then again, I don't expect to see much noise in that range and don't want to needlessly distort the signal. Should I just send the latest sample, or perhaps take the average of all samples between two "PC sample points"? What would you recommend?

Edit to address Oscar's questions:

I don't think there will be much signal above 5Hz and the SNR should be pretty good without any filtering already, but I'll try recording some examples soon. Different environments and different source equipment for the signal could also look quite different, though.

At the moment, a different device that we use for the same purpose samples the input at only 2Hz without extra filtering after the A/D, but the A/D module in that device seems to filter quite a bit already. Trying to get the requirements, I got the "at most 5Hz" figure from a colleague (though I'm not even sure if he meant samplerate or nyquist frequency), so I adjusted the sample rate up in order to remain consistent.

I'll probably have to handle four channels sampled at 12 bits, and the filtering is the only CPU-intensive task the controller has to perform. The number of cycles per sample is a tradeoff against sample rate (ammount of aliasing) and code maintainability (perform calculations while communicating with the ADC?), but I think 400 cycles would be a pretty good result.

I can't push the full sample rate to the PC, but 100-200sps per channel should be ok.

Is this a good filter for my task, or would you recommend something else?

It might be given that you do not need to get rid of any other noise/signal in that band. The advantage will be the lack of multiplications and therefore most likely fewer cycles in your implementation.

  Should I just send the latest sample, or perhaps take the average of all samples between two "PC sample points"?

Taking the average is just like using another moving average filter, so this may or may not be a good idea. It is easily analyzed in a similar way. Given that your signal is "good enough" to downsample, filtering again will, while removing a bit of noise, primarily mess up with the passband.

It is worth noting though that it may be worthwhile to use a so called polyphase decomposition. You do not really want to spend time computing intermediate samples which are not transferred to the PC anyway. Hence, you do in fact have quite a bit more cycles per sample. The straightforward way to go about this to use an (or several) FIR filter. In the case of several compute the convolution of all. Only apply the filter such that it computes exactly the output you want. If you go with this approach you might just as well design a single FIR filter using general multipliers (since the convolution of the moving average filters will have general coefficients anyway).

Assuming you have the signal package for Octave (and that the call is identical to Matlab), you can try something like:

fs=1000
fn = fs/2 % To simplify things
delta = 3 % range in Hz around the power line frequencies, should be less than 5...
eps = 1e-7 % Some small number....
N = 50 % Filter order
W = 10 % How much more you'd like to attenuate the power compared to generic noise

h = firpm(N, [0 5 (50-delta) (50+delta) (50+delta+eps) (60-delta-eps)  (60-delta)  (60+delta) (60+delta+eps) fn]/fn,[1 1 0 0 0 0 0 0 0 0], [1 W 1 W 1])

The lesser you pick delta, the more likely that you will get a notch exactly where you want it...


Line Detection with OpenCV and Python Why do we get less amount of lines then desired in line detection using HoughLines tranformation from OpenCV? How could we improve this method and get all the lines appearing in the image? I would still use HoughLines if probable.

The performance of line detectors are image and application dependent in generally. Typically, people tune the parameters of the detectors to get the best possible result. When such a result cannot be obtained, it is often because of the insufficient performance of the detector for the current problem. One approach would be to change the algorithm. For instance, OpenCV also has an LSD line segment detector. Have you tried that? Maybe it can also be used in combination with Hough lines.

Having said all those, I would also look at the hardware/system related issues, if they are under your control: Smartly selecting system parameters such as lighting, cameras and lenses generally solve half of the problem. If you would like to operate in the wild, then I would suggest more sophisticated software solutions.

Edit: This question looks very similar to the one recently posted:
MatLab - Can not detect all lines in a "simple" image using houghlines

Maybe this could help you out.


What is the relationship between the sigma in the Laplacian of Gaussian and the two sigmas in the Difference of Gaussians? I understand that a Laplacian-of-Gaussian filter can be approximated by a Difference-of-Gaussians filter, and that the ratio of the two sigmas for the latter should be 1:1.6 for the best approximation. However, I'm not sure how the two sigmas in the Difference of Gaussians relates to the sigma for the Laplacian of Gaussian. Is the smaller sigma in the former equal to the sigma of the latter? Is the larger sigma? Or is the relationship something else?

I understand that a Laplacian-of-Gaussian filter can be approximated by a Difference-of-Gaussians filter, and that the ratio of the two sigmas for the latter should be 1:1.6 for the best approximation

In theory, the smaller the ratio between two sigmas, the better the approximation. In practice, you'll get numerical errors at some point, but as long as you're using floating point numbers, smaller values than 1.6 will give you a better approximation.

To illustrate, I've plotted a cross-section of the LoG and DoG for a few values of k in Mathematica:

As you can see, k=1.6 is not an ideal approximation. For example, k=1.1 would give a much closer approximation.

But you usually want to calculate LoG approximations for a range of sigmas. (Otherwise, why bother with the DoG approximation at all? Calculating a single LoG filtered image isn't more expensive than calculating a single DoG filtered image.) So the value of k is usually chosen so that you can calculate a series of gaussian filtered images with sigmas s, sk, sk^2, s*k^3..., and then calculate the differences between adjacent gaussians. So if you choose a smaller k, you'd have to calculate more "layers" of gaussians for the same sigma-range. k=1.6 is a trade-off between wanting a close approximation and not wanting to calculate too many different gaussians.

  However, I'm not sure how the two sigmas in the Difference of Gaussians relates to the sigma for the Laplacian of Gaussian. Is the smaller sigma in the former equal to the sigma of the latter?

From the formulas on the wiki page @Libor linked to, you can see that $t=\sigma ^2$, so the approximate a LoG for some sigma, you need two gaussians with sigmas $\sqrt{\sigma ^2+\text{$\Delta $t}}$ and $\sqrt{\sigma ^2-\text{$\Delta $t}}$ (at least in the limit $\text{$\Delta $t}\to 0$). Or, in terms of k:

$\sigma _{\text{Laplace}}=\sigma \sqrt{\frac{1+k^2}{2}}$


FFT: low magnitude peaks centered exactly 1 Hz apart across whole spectrum I am running an fft on 60 seconds of data sampled at 1289Hz. The fft is run every second on a single channels fifo buffer which is also updated every second. The algorithm used is numpy rfft with a numpy blackman window. The fft output is normalized to the max fft value. This is all on a raspberry pi 4 with custom adc board that uses an AD4114. Am sampling 8 channels (16 in differential mode) as fast as the AD4114 can go which works out to 10.312kHz.
The following screen shots are with a 20Hz sign wave applied to one of the channels.  Similar low magnitude peaks appear with real data a transducer attached to the input. There are low magnitude peaks centered exactly 1Hz apart right on top of each whole number. Is this normal? Like a normal artifact of the fft? Possibly something wrong with my data? I am a casual fft user and am hoping someone may have some advice for me. Thank you very much : )

Is this normal?

No

Like a normal artifact of the fft?

No

Possibly something wrong with my data?

Yes. Could be a data problem, could be a hardware problem, could be a SW problem, could be real. You seem to have a very large DC offset in your data. Without knowing what exactly your data, is I can't tell whether that's expected or if it's a problem.


Skeletonization of a Circle - Image Processing What will be the output of skeletonization of a circle?
In my class my instructor said that it will be a single dot in the center but what I found on the Internet was a cross.
I discussed with a friend of mine, and he said that the reason it will be a dot is that a smaller circle can not fit a larger circle in such a way that it touches more than one point and is not out of the boundary of the original signal.
Although I beg to differ because of the equation of finding the length of the arc.
Can anyone confirm the output of skeletonization of a circle with proper reasoning?

The question here is: "The output of which skeletonisation"? The topological or the morphological?

It seems from what you mention on your discussion with your friend that you are referring to the morphological skeleton.

In the case of the topological skeleton, the answer is trivial. Given that the topological skeleton is the set of all foreground pixels that are at maximal distances from the boundaries of the image and the definition of the disk (all points in space that are at most at distance $R$ from the centre of the disk), the skeleton of a disk (a circle that is filled*) is the single pixel at its centre.

The morphological skeleton works slightly differently. It is based on the morphological operations of erosion followed by opening, both of which are using a structuring element.

Roughly speaking, you can think of erosion followed by opening as progressively "pealing off" boundaries from a region, until there is nothing left to peal off. (This is not entirely accurate and I will return to it further below).

Here for example is some random disc I placed on a dark background...

...and here are a few of the first iterations of morphological skeletonisation:

Where it is clear that the disc is shrinking and this looks like going towards our previous thinking regarding the topological skeleton, where eventually, all that is left is going to be a single pixel at the centre.

However, we can already see that this shape does not look so much like a circle now, it's kind of pointy and also, if we continue this process, something odd starts happening (in my case after about iteration 12 or 13):

Ultimately, what is left behind is a cross.

Where is this cross coming from?

This is a byproduct of two things:

The fact that the structuring element is discrete and its curvature at different iteration levels does not match exactly the curvature of the given disc. To understand this better, consider this extreme example that could well be showing a "disc" or a "diamond" or a square rotated 45 degrees:

$$\begin{matrix}
    0 & 0 & 1 & 0 & 0 \\
    0 & 1 & 1 & 1 & 0\\
    1 & 1 & 1 & 1 & 1\\
    0 & 1 & 1 & 1 & 0\\
    0 & 0 & 1 & 0 & 0\\
\end{matrix}$$

(In other words, the disc is not exactly a disc and it is not strictly enforcing that rule of "all pixels at distance $R$" strictly to all of its pixels. It is more like, pixels less than that distance in which case, its discrete boundary "plays" in and out of the actual continuous boundary.)

The way that the "pealing off" takes place. In fact, the morphological skeleton is the union of $n$ partial skeletons. Each partial skeleton of some image $I$ via a structuring element $B$ is obtained via:

$$ S_n(I) = (I \ominus nB) - (I \ominus nB) \circ B $$

Where $I$ is the image (here, our bigger disc), $n$ is the iteration level, $B$ is the structuring element, $\ominus$ represents erosion, $\circ$ represents dilation and $nB$ is interpreted as "$n$ applications of erosion". Also notice here that the $-$ operator has the lowest precedence. So you don't evaluate it until you have evaluated both the left-hand and right-hand expressions.

The "pealing off" is achieved by that $-$ operator in the middle but there is a little bit more into it. A partial skeleton is the difference (as in  set difference) between the "eroded" $I$mage by $B$ and the "dilation of that eroded $I$mage" by $B$. Which is basically like the difference between a shrunk image and the shrunk image enlarged. This difference is $\emptyset$ for the first few iterations but as the disc gets smaller and smaller and now the differences in the boundary start becoming more pronounced, the shrunk and shrunk-enlarged images start differing by 1 pixel (exactly at the peaks as you can see from the matrix example above). Since the skeleton is the union of the partial skeletons and each partial skeleton will have 4 pixels difference to its previous one, what "survives" throughout the iteration levels is that cross.

To understand better why I am referring to erosion and dilation as shrinking and enlarging and their difference as some "boundary", you need to understand the effect of those particular operations in detail. To do this, in addition to the excellent wikipedia articles, you can have a look at Gonzalez & Wood's "Digital Image Processing". In my second edition, chapter 9 is devoted to morphological image processing (and is totally worth a complete thorough reading) and the morphological skeletons are presented in pages 543-545.

Hope this helps.

*: The topological skeleton of a circle is a circle.


What is the advantage of MATLAB's filtfilt MATLAB's filtfilt does a forward-backward filtering, i.e., filter, reverse the signal, filter again and then reverse again. Apparently this done to reduce phase lags? What are the advantages/disadvantages of using such a filtering (I guess it would result in an effective increase in filter order).

Would it be preferable to use filtfilt always instead of filter (i.e., only forward filtering)? Are there any applications where it is necessary to use this and where it shouldn't be used?

You can best look at it in the frequency domain. If $x[n]$ is the input sequence and $h[n]$ is the filter's impulse response, then the result of the first filter pass is

$$X(e^{j\omega})H(e^{j\omega})$$

with $X(e^{j\omega})$ and $H(e^{j\omega})$ the Fourier transforms of $x[n]$ and $h[n]$, respectively. Time reversal corresponds to replacing $\omega$ by $-\omega$ in the frequency domain, so after time-reversal we get

$$X(e^{-j\omega})H(e^{-j\omega})$$

The second filter pass corresponds to another multiplication with $H(e^{j\omega})$:

$$X(e^{-j\omega})H(e^{j\omega})H(e^{-j\omega})$$

which after time-reversal finally gives for the spectrum of the output signal

$$Y(e^{j\omega})=X(e^{j\omega})H(e^{j\omega})H(e^{-j\omega})=
X(e^{j\omega})|H(e^{j\omega})|^2\tag{1}$$

because for real-valued filter coefficients we have $H(e^{-j\omega})=H^{*}(e^{j\omega})$. Equation (1) shows that the output spectrum is obtained by filtering with a filter with frequency response $|H(e^{j\omega})|^2$, which is purely real-valued, i.e. its phase is zero and consequently there are no phase distortions.

This is the theory. In real-time processing there is of course quite a large delay because time-reversal only works if you allow a latency corresponding to the length of the input block. But this does not change the fact that there are no phase distortions, it's just an additional delay of the output data. For FIR filtering, this approach is not especially useful because you might as well define a new filter $\hat{h}[n]=h[n]*h[-n]$ and get the same result with ordinary filtering. It is more interesting to use this method with IIR filters, because they cannot have zero-phase (or linear phase, i.e. a pure delay).

In sum:

if you have or need an IIR filter and you want zero phase distortion, AND processing delay is no problem then this method is useful
if processing delay is an issue you shouldn't use it
if you have an FIR filter, you can easily compute a new FIR filter response which is equivalent to using this method. Note that with FIR filters an exactly linear phase can always be realized.


Changing camera and image properties to calibrate a camera Is it OK to calibrate a camera using effects and filters such as increasing contrast and brightness etc. to improve the detection of corners of a chessboard?
Even if the used mode for detecting the chessboard is just plain gray scale?

You should not experience any major drawbacks.
Generally such techniques utilize adaptive thresholds and rely on gradient orientations (especially normalized ones), which are intrinsically invariant to uniform illumination and contrast changes. In other words, brightness changes will not effect your subpixel positions, if you are using an appropriate corner detection algorithm.


Finding number of alternations in parks mcclellan I'm given this question and asked to find number of alternations for filter 1 and 2, using parks mcclellan algorithm.

I search for solution and find this answer :

My question is why the number of filter taps are 9 in the solution, and how exactly to use parks mcclellan to resolve number of alternations for filter 1 and 2 ?

This question is from Oppenheim, Schafer and Buck's textbook Discrete-time Signal Processing. There must be a typo either in the problem or in the solution manual. But the only important thing is to be able to find the number of alternations, and to be able to check whether a given filter can be optimal for a given filter length.
For a type I (even symmetry, odd filter length) linear phase low pass filter, the degree of the corresponding polynomial is $L=(N-1)/2$, where $N$ is the filter length. Note that the filter order $M$ is equal to $N-1$.
According to the alternation theorem, a Chebyshev optimal type I low pass filter must have $L+2$ or $L+3$ alternations. Alternations are points of maximum error, with the only restriction that adjacent alternations must have errors of different signs (that's why we talk about alternations after all). Note that the band edges also count as alternations because at these points the approximation error attains its maximum.
Looking at the figures we see that the first filter has $3$ alternations in the passband (one at $\omega=0$, one at the maximum, and one at the passband edge), and $3$ alternations in the stopband (one at the band edge, one at the minimum, and one at $\omega=\pi$), which makes for a total of $6$ alternations. In the same way you can see that the second filter has $7$ alternations.
If the filter length is $N=11$, we get $L=(N-1)/2=5$, and the consequently, the optimal filter must have $L+2=7$ or $L+3=8$ alternations. In this case, only the second filter satisfies the alternation theorem.
On the other hand, for $N=9$ we have $L=4$, and both filter satisfy the alternation theorem because $L+2=6$ and $L+3=7$.


Medical image dataset I am doing some project on medical image processing and I need some uncompressed medical images especially magnetic resonance angiography, vessel and so on. I have searched a lot but most of the available images are in JPEG format which is not a suitable format due to its lossy compression. Do you know any free dataset which have uncompressed medical images especially those mentioned?

The problem you are most likely having is that you need to search for files in the correct format.  DICOM is what you need to look for.
A few hits:

One example
A library from a DICOM viewer software producer
A library of images from Siemens

You will also need a DICOM viewer.
Image J  is free (non-commercial, no restrictions on use) and can be used under Windows, Mac, and Linux - it does require Java.
There are other viewers available, free as well as commercial.  Just search for DICOM viewers.


Confused in difference between sequence versus signal? I am learning basics of signals and systems with MATLAB
If i create a vector in matlab such that
x=[ 2 5 10 20], will it be considered a sequence or a signal?

As far my learning/knowledge is, x is a sequence in this case, since it is using discrete numbers as its elements.

If there was x such that x=sin(t) or cos(t) or exp(t),it would be a signal because of the continuous range of values in it.

Please correct/guide me if i am wrong

In DSP we use the term signal to denote continuous-time functions as well as discrete sequences. To avoid confusion you can use the terms continuous-time signals and discrete-time signals.

Of course, the independent variable doesn't have to be "time" (e.g., in image processing). In that case you would simply be talking about continuous or discrete signals.

Note that often you would come across the terms analog or digital signals for continuous-time and discrete-time. That categorization is a bit more ambiguous because by digital we usually mean that not only the independent variable (often "time") is discrete but also the signal values, i.e., that the signal is quantized. But in practice it's not always clear if digital is used in that strict sense. So if you don't want to restrict yourself to quantized sequences, it is advisable to use discrete(-time) instead of digital.

Take a look at this question and its answers for more details on the distinction continuous/discrete and analog/digital.


Very simple problem regarding shifting of a signal I have this simple $x(t)$ graph, and I have some confusion regarding time shifting.

What is the graph of $x(1-t)$ ?

Don't look at graphs and don't think of "flip and slide" as your
instructor/TA will tell you is the way to solve such problems.

$x(\cdot)$ is a function that has value $5$ whenever what
is between the "$($" and the "$)$" has value between $-4$ and $4$.  So, $x(1-t)$
has value $5$ exactly when

$$ -4 \leq 1-t \leq 4 \quad \Rightarrow \quad t \in [a, b]$$

where $a$ and $b$ are numbers (hint: they are integers!) that I will
leave for you to work as an exercise. When you are done, you can easily
sketch $x(1-t)$  from this information.


While finding the ESD/PSD of a signal why we always prefer to find it via Auto-correlation function then the square of the FT of the signal? In a video i saw that while calculating ESD or PSD of a signal time auto correlation function was used when it can be also done by the square of FT of the signal.Why we followed that approach even when the later process is more straight-forward and easy.

How does it helps us in later stage?

In the theory of discrete-time random signals, the spectral description of a random process $X[n,s]$ (abbreviated $x[n]$) is often desired, just as it's desired for deterministic signals.
When the frequency of a signal is mentioned at once, the immediate consequent is the Fourier transform of that signal. However, when the signal is random, there is a problem: random processes do not have Fourier transforms. They have infinite energy as the sum
$$ \sum_{n=-\infty}^{\infty} |x[n]|^2 $$
diverges; (as $\lim_{n \to \infty} |x[n]|^2 \neq 0$, which is sufficient to prove that the infinite series cannot converge.) . Thus, the DT Fourier transform of $x[n]$ does not exist:
$$X(\omega) = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n} $$ does not converge for any $\omega$.
To overcome this problem for wide-sense stationary random processes, we define and use the auto-correlation sequence (ACS) $r_{xx}[m]$ of the (zero-mean) random process:
$$ r_{xx}[m] = \mathcal{E}\{x[n]x[n+m]^{*} \} ~~~,~~~\text{ for } -\infty < m < \infty$$
The most important property of this ACS (from DTFT point of view) is that it has finite energy (assume aperiodic case) which is sufficient to compute a Fourier transform of it denoted as the power spectral density (PSD) of the WSS random process $x[n]$:
$$ S_{xx}(\omega) = \sum_{m=-\infty}^{\infty} r_{xx}[m] e^{-j \omega m} $$
All these are theory. Because to compute the PSD, you need to know the ACS. But the ACS is related to the PDF (probability density/distribution function) of the random process $x[n]$. And the PDF in practice is not known. In theory (as a part of a homework assignment for example) you can propose a PDF for a given process, and then proceed analytically to compute its ACS and PSD. Yet in practice the signals that you are observing are truly random hence the analytical formulation of their PDF is not exactly known.
So in practice you must find a convenient means of computing the PSD of a partially observed random process (such as a speech given my Martin Luther King). All these approaches are rooted in what is called as Estimation Theory which deals with finding convenient methods of computing certain statistical parameters related with random variables and processes.
One of the methods to estimate the PSD of a random process, relies on estimating its ACS from the observed data, and then taking the Fourier transform of it to reach an estimate of the true PSD. (The true power spectral density relies on the true ACS which is not available but estimated from data).
The other approach, and the more commonly used one, directly estimates the PSD (which is itself a random quantity) from the observed data. The most basic of such direct PSD estimation is known as the periodogram and obtained by computing the Fourier transform of the signal itself;i.e.,
$$I(\omega) = K \cdot |X(\omega)|^2 $$ where $X(\omega)$ is the DTFT of the finite length and finite energy observed random signal $x[n]$,  $K$ is a proper scaling factor and $I(\omega)$ is the periodogram estimate of the PSD.
And this is shows a practical means of computing (actually estimating) the PSD associated with a random process $X[n,s]$ over its partial observations $x[n]$.


Filter Order Rule of Thumb Is there a rule of thumb for calculating the approximate filter order of a FIR/IIR based on passband vs sample rate?  I seem to remember reading about it once in Undergrad. I'd think you could make a rough estimate like scaler * passband/stopband.

I have been working on implementing different FIR/IIR filters and would like some rough formula that says if one variable is changed, this is approximately how many more taps you would need to get similar performance.

My favorite "Rule of thumb" for the order of a low-pass FIR filter is the "fred harris rule of thumb":
$$
N=\frac{f_s}{\Delta f}\cdot\frac{\rm atten_{dB}}{22}
$$
where

$\Delta f$ is the transition band, in same units of $f_s$
$f_s$ is the sample rate of the filter
$\rm atten_{dB}$ is the target rejection in dB

For example if you have a transition band of 100 Hz in a system sampled at 1 kHz, and your rejection requirement is 50 dB in the stop band then the order can be approximated by:
$$
N = \frac{1\ \rm kHz}{100\ \rm Hz} \cdot\frac{50}{22} = 23\ \rm taps\tag{rounding up}
$$
Thank you Fred Harris!
Note another more detailed formula that takes into account passband ripple is Kaiser's formula thanks to James Kaiser from Bell Labs, which I included in my graphic below.
For most applications I have done, the Fred Harris approach has been just fine, as given a certain rejection, the resulting filters using traditional filter design algorithms such as Parks-McClellan and Remez have exceeded my passband ripple requirements upon meeting the rejection requirement.  (What I usually do is estimate the order, design the filter with that order, inspect the result and increase or decrease order from there to fine tune). The results from the estimates are just that: estimates, and can vary widely depending on overall design parameters and not to be assumed to be an exact solution.

For those familiar with filter design using window approaches, review of the boxcar or rectangular window (which is simple truncation) reveals why it takes $f_s/\Delta f$ taps (which is the same as $2\pi/\Delta\omega$ if the units for normalized frequency are radians/sample as is often done) to complete the transition band. See the images below that help explain this.
The top image below shows the expected Sinc in frequency for a rectangular window in time, in this case as a non-causal rectangular pulse centered on $t=0$. This is then repeated in discrete forms as a causal waveform starting at $t=0$, both with the Discrete Time Fourier Transform (DTFT) and Discrete Fourier Transform (DFT) where the difference is the samples in time extend to $\pm \infty$ for the DTFT resulting in a continuous waveform in the frequency domain. In both cases the result is an aliased Sinc function that is periodic over the interval $f= [0,f_s)$, with the key point that for $N$ samples in time of the rectangular function, the frequency response will have its first null at $f = 1/N$ (Where $f$ is the normalized frequency with 1 being the sampling rate).

This next image below shows the rectangular window approach to filter design (which I would never recommend but is informative). The first plot in the upper left hand corner shows the target frequency response for our filter as an ideal "brick wall" response. Please don't confuse that with the "boxcar window" (or "rectangular window") which is also a rectangular shape-- the window is in the time domain!
In order to realize such a filter, we would use the impulse response of the desired frequency response as the coefficients in our FIR filter (the coefficients of the filter are the impulse response---put an impulse in and out come all the coefficients!).  The impulse response for a rectangular frequency (brickwall) response is the inverse FT which is a Sinc function, in the time domain, shown in the lower left hand corner as the "Required Impulse Response". A Sinc function extends to plus and minus infinity, so to actually realize such a filter we would need an infinitely long filter and it would have infinitely long delay. Obviously we can't do that, so we truncate the coefficients to something realizable. The longer the filter, the closer we approximate the ideal brickwall response, but also the longer the delay will be (and the more resources we would need in terms of the filter construction; more taps).
Truncating the impulse response in the time domain is mathematically identical to multiplying by a rectangular window in the time domain. (Note that the impulse response is also delayed by half the duration of the window in order for the system to be causal).  Multiplying in the time domain is equivalent to convolution in the frequency domain. The frequency domain (FT) of the impulse response prior to truncation is our original desired brickwall frequency response. The frequency response for the rectangular window is a Sinc function in the frequency domain.
So when we truncate the desired impulse response (multiply in time by a rectangular window), we convolve the desired frequency response with a Sinc function, resulting in an approximation of our target frequency response as shown in the upper right-hand corner of the image below.

A key take-away for Sinc functions in general is the first null is $1/T$ where $T$ is the duration of the rectangular function.  For a sampled system, the first null would be at $2\pi/N$ where $N$ represents the number of samples for the duration of the rectangular function.  In the images a normalized radian frequency is used for the frequency axis- (if that confuses you just know the $2 \pi$ is the radian frequency for the sampling rate).  So in the process of the convolution, the sharp brickwall transition spreads out and in this case goes to 0 (our $\Delta \omega$) over a frequency of $2\pi/N$!  So here $$N = 2\pi/\Delta \omega$$ and of course the filter is poor with sidelobes etc.  Note this: This transition from the Sinc function is the sharpest available for a given number of taps; it has the best resolution in frequency but the poorest dynamic range (rejection). Other window typologies (Blackman, Blackman-harris, Kaiser (my favorite), etc) will significantly improve dynamic range but always at the expense of the transition.

So from the above, we see the origin of the $2\pi/\Delta\omega$ that is used in the approximation formulas, and we also see why there is an additional multiplication factor increasing the number of taps above this for typical filter designs; the rectangular window would give us the best transition possible with $N$ taps where $N = 2\pi/\Delta\omega$ but have very poor rejection. More taps are used to smooth the time transition further beyond the sharp transition of the rectangular window, providing greater rejection at the expense of transition bandwidth.


Digital filter design basic principles (IIR/FIR) Although I have a solid experience in designing audio engines and such, I am fairly new to the realm of Digital Filter Design, particularly IIR and FIR filters. In other words, I'm trying to learn as much as I can on how to design filters and derive their difference equations. I'm starting from the basics, so please bear with me, like I said, I'm trying to learn.

Herein is my question:

Say I want to design a low-pass filter with a particular cutoff - say 300 Hz. What would be the best way of mathematically deriving the transfer function and then deriving the particular difference equation to implement the filter in Direct Form I and Direct Form II (or only DF-1 for now...)?

I have some understanding of transfer functions and how they relate to difference equations from some excellent material on the web, unfortunately some of it assumes a good bit of prior knowledge, so it is more confusing than helpful in my quest. So I guess I need a more step-by-step example which will help me connect the dots.

So I'm basically looking for help with a breakdown of the process starting from choosing the cutoff frequency up to deriving the difference equation.

Any help will be much appreciated. I am familiar with a lot of the concepts - impulse response, DFT's, the math behind it, I guess what I need more help with is the concept of the z-transform and the pole/zero use to design the filter's transfer function and then how does the cutoff freq. play in all this to finally derive the difference equation.

Since I tend to learn best from example, I thought I'd ask here.
Thanks a lot to anyone who finds the time to help me out.

Digital filter design is a very large and mature topic and - as you've mentioned in your question - there is a lot of material available. What I want to try here is to get you started and to make the existing material more accessible. Instead of digital filters I should actually be talking about discrete-time filters because I will not consider coefficient and signal quantization here. You know already about FIR and IIR filters, and you also know some filter structures like DF I and II. Nevertheless, let me start with some basics:

A non-recursive linear time-invariant (LTI) filter can be described by the following difference equation

$$y(n)=h_0x(n)+h_1x(n-1)+\ldots +h_{N-1}x(n-N+1)=\sum_{k=0}^{N-1}h_kx(n-k)\tag{1}$$

where $y(n)$ is the output sequence, $x(n)$ is the input sequence, $n$ is the time index, $h_k$ are the filter coefficients, and $N$ is the filter length (the number of taps). The filter taps $h_k$ are also the impulse response of the filter because if the input signal is an impulse, i.e. $x(n)=\delta(n)$, then $y(n)=h_n$ (if the filter's memory has been initialized with zeros). Equation (1) describes a linear time-invariant finite impulse response (FIR) system. The sum on the right-hand side of (1) is a convolution sum, i.e. the output signal is obtained by convolving the input signal with the impulse response. This is always true, but for IIR filters we cannot explicitly compute the convolution sum because the impulse response is infinitely long, i.e. there are infinitely many coefficients $h_k$. One important advantage of FIR filters is that they are always stable, i.e. for a bounded input sequence, the output sequence is always bounded. Another advantage is that FIR filters can always be realized with an exactly linear phase, i.e. they will not add any phase distortion apart from a pure delay. Furthermore, the design problem is usually easier, as we will see later.

A recursive LTI filter is described by the following difference equation:

$$y(n)=b_0x(n)+b_1x(n-1)+\ldots+b_Mx(n-M)-\\
-a_1y(n-1)-\ldots-a_Ny(n-N)\tag{2}$$

Equation (2) shows that the output is not only composed of weighted and delayed input samples, but also of weighted past output samples. In general, the impulse response of such a system is infinitely long, i.e. the corresponding system is an IIR system. However, there are special cases of recursive filters with a finite impulse response. Note that the impulse response is not anymore given by either the coefficients $b_k$ or $a_k$ as in the case of FIR filters. One advantage of IIR filters is that steep filters with high stopband attenuation can be realized with much fewer coefficients (and delays) than in the FIR case, i.e. they are computationally more efficient. However, one needs to be careful with the choice of the coefficients $a_k$ because IIR filter can be unstable, i.e. their output sequence can be unbounded, even with a bounded input sequence.

Filters can be designed according to specifications either in the time (sample) domain or in the frequency domain, or both. Since you've mentioned a cut-off frequency in your question, I assume you're more interested in specifications in the frequency domain. In this case you need to have a look at the frequency responses of FIR and IIR systems. The frequency response of a system is the Fourier transform of its impulse response, assuming that it exists (which is the case for stable systems). The frequency response of an FIR filter is

$$H(e^{j\theta})=\sum_{k=0}^{N-1}h_ke^{-jk\theta}\tag{3}$$

where $\theta$ is the discrete-time frequency variable:

$$\theta=\frac{2\pi f}{f_s}$$

with the actual frequency $f$ and the sampling frequency $f_s$. From (3) you can see that approximating a desired frequency response by an FIR system is basically a problem of polynomial approximation. For recursive systems we have

$$H(e^{j\theta})=\frac{\sum_{k=0}^Mb_ke^{-j\theta}}{1+\sum_{k=1}^Na_ke^{-j\theta}}\tag{4}$$

and you get a rational approximation problem, which is usually much more difficult than the polynomial approximation problem in the case of FIR filters. From (3) and (4) you can see that the frequency response of an FIR filter is of course only a special case of the response of a recursive filter with coefficients $a_k=0$, $k=1,\dots,N$.

Let's now take a quick look at filter design methods. For FIR filters you could take an inverse Fourier transform of the desired frequency response to get the impulse response of the filter, which directly corresponds to the filter coefficients. Since you approximate the desired response by a finite length impulse response you should apply a smooth window to the obtained impulse response to minimize oscillations in the actual frequency response due to Gibbs' phenomenon. This method is called frequency-sampling method.

For simple standard filters like ideal lowpass, highpass, bandpass or bandstop filters (and a few others), you could even analytically calculate the exact impulse response by taking the inverse Fourier transform of the ideal desired response:

$$h_k=\frac{1}{2\pi}\int_{-\pi}^{\pi}H(e^{j\theta})e^{jk\theta}d\theta$$

This integral is easy to evaluate for piecewise constant desired responses, as is the case for ideal frequency-selective filters. This will give you an infinitely long, non-causal impulse response, which needs to be windowed and shifted to make it finite and causal. This method is know as window-design.

There are of course  many other FIR filter design methods. One important numerical method is the famous Parks-McClellan exchange algorithm which designs optimal filters with constant passband and stopband ripples. It is a numerical approximation method and there are many software implementations available, e.g. in Matlab and Octave.

The most common IIR design method for frequency selective filters is the bilinear transformation method. This method simply uses analytical formulas for the design of optimal analog filters (such as Butterworth, Chebyshev, Cauer/elliptic, and Bessel filters), and transforms them to the discrete-time domain by applying a bilinear transformation to the complex variable $s$ (analog domain) which maps the (imaginary) frequency axis of the complex $s$-plane to the unit circle in the complex $z$-plane (discrete-time domain). Don't worry if you do not yet know much about complex transfer functions in the analog or discrete-time domain because there are good implementations available of the bilinear transform method, e.g. in Matlab or Octave.

There are of course many more interesting and useful methods, depending on the type of specifications you have, but I hope that this will get you started and will make any material you come across more understandable. A very good (and free) book covering some basic filter design methods (and a lot more) is Intoduction to Signal Processing by Orfanidis. You can find several design examples there. Another great classic book is Digital Filter Design by Parks and Burrus.


What is difference between Signal frequency and data rate I read somewhere the below two line and got confused what is the actual difference between the two line. These two lines are related to CAN

1.Maximum Signal frequency used is 1 Mbit/sec

2.Typical maximum data rate achievable is 40KBytes/sec.

If anyone knows about it please share the answer.

There is a mistake. The units of frequency is cycles/second and not Mbit/second. I noticed this mistake in one of the online links on CAN.

I think what it means is  that the maximum theoretically possible data rate is 1Mbps and the typical real life maximum data rate is 320Kbps(40KBytes/sec). The standard says that for network lengths below 40m, 1Mbps can be achieved but for network lengths of 500m the data rate achieved is around 125kbps. So I think that's reason typical data rate is given as 40KBytes/sec.

References:

http://en.wikipedia.org/wiki/CAN_bus
http://www.computer-solutions.co.uk/info/Embedded_tutorials/can_tutorial.htm


Frequency Spectrum Analysis How do you find the phase and inverse Fourier transform of this frequency spectrum?

The magnitude is given by $ {\left| X \left( f \right) \right|} $ hence it is 6.
Assuming all data is given in the graph, the function isn't complex it is only negative -> Phase is -180 degrees.
You should use its symmetry. Without over thinking it, It could be generated by a window on the center convolved with 2 deltas (Shifted) -> Sinc in time multiplied by 2 exponentials.


Magnitude of Sobel filter result, vector or scalar? I have a question about image gradient magnitude, given an image $I$. The horizontal and vertical derivatives resulting from Sobel filters are defined as:
$$K_x=
        \begin{bmatrix}
        -1 & 0 & 1 \\
        -2 & 0 & 2 \\
        -1 & 0 & 1 \\
        \end{bmatrix}*I;
K_y=
        \begin{bmatrix}
        -1 & -2 & -1 \\
         0 & 0 & 0 \\
         1 & 2 & 1 \\
        \end{bmatrix}*I
$$

where $*$ denotes the convolution operator. The gradient magnitudes can be computed as

$$I_{\delta}=\sqrt{(K_x^2+K_y^2)}$$.

Now, I want to "invert" the result of $I_{\delta}$, so that maximum magnitude gradient values become zero. Then the "inverting" function $g$ is
$$g(I_{\delta})=\max (I_{\delta}) - I_{\delta}\,.$$

My question is: is $I_{\delta}$ scalar or vector?
What do you think? If it is possible, could you correct help me.

Pretty much by definition, the magnitude of anything is a scalar value.  Though, since you are working on an image, the result is a 2D array of scalar values.

Given Kx and Ky, it's also possible to work out the directions of the edges, as well as the magnitudes, but that's a whole different bit of maths.


Cross correlation and cross power spectrum for signal + noise Quick question about my understanding of cross correlations and cross spectral densities.

Let $C_{ab}$ be the cross correlation between the signals $a$ and $b$, i.e., $C_{ab}(\tau) = \langle a(t)b(t+\tau)\rangle$. Let $S_{ab}$ be the cross spectral density of the signals $a$ and $b$, i.e., $S_{ab}=\mathcal{F}(a)^* \mathcal{F}(b)$ and by the convolution theorem/cross correlation theorem $S_{ab} = \mathcal{F}(C_{ab}(\tau))$.

Say $x(t)$ is some signal, and $\nu(t)$ is some uncorrelated noise signal. Let $s(t)=x(t)+\nu(t)$
Thus, $C_{xs} = C_{xx}+C_{x\nu}$ from the linearity of the expectation value being calculated in the cross correlation. By initial assumption that $x$ and $\nu$ are uncorrelated, we get $C_{xs}=C_{xx}$. (This calculation so far also appears here on page 5, for example so I think it is correct).

Now, if we take the fourier transform of both sides of this, we get $S_{xs}=S_{xx}$. Is this correct so far? But does this not mean that $\mathcal{F}(x)^* \mathcal{F}(s) = \mathcal{F}(x)^*\mathcal{F}(x)$? and hence $\mathcal{F}(s)=\mathcal{F}(x)$, and hence even $s(t)=x(t)$? Clearly there is something wrong here. What mistake am I making?

Thanks!

EDIT: On a similar note, if $C_{x\nu}=0$, does that not mean that $\mathcal{F}(x)^*\mathcal{F}(\nu)=$ and hence (by taking the absolute value of both sides and squaring) that $S_{xx} S_{\nu\nu}=0$? This is patently false. At each frequency this would require either the power in my signal to be zero, or the power in the noise to be zero. Even if my noise is white noise then this is clearly false. I'm doing something fundamentally wrong here. What? Any help?

$S_{xy}=S_{xx}$ is true but this does not mean $F_x*F_y = F_x*F_x$ when $y$ is random.

The actual definition is $$S_{xy}(f)=\lim_{T\rightarrow \infty}\textbf{E}\{F^*_xF_y\}$$ for random signals. Since $x$ and $v$ are uncorrelated, this would be equal to $$\lim_{T\rightarrow \infty}\textbf{E}\{F^*_xF_x\} = S_{xx}.$$ But your next assumption is not correct.

$F_x$, $F_y$ are only based on single realization of $x$ and $y$ respectively. For random signal, $S_{xy} \ne F_x^*F_y$


normalized frequency to frequency as the definition of normalized frequency is the frequency divided by the sample frequency, I'm confused about my signal so as an example a signal sampled at ($f_s=500 000$) Hz when applying a FFT, an impulse at frequency of $87$ Hz occurred(true value confirmed) while it occurred at $0.117$ Hz in normalized frequency so we know that 1 for normalized frequency is sampling frequency($f_s$) but $0.117$ Hz normalized frequency how it calculated according to $f_s$ to get $87$ Hz while it's not $87/(f_s/2)$!!
I will be grateful for your help

You have to carefully read the documentation of the functions that use/generate normalized frequency in some cases they use $fs$ and in others they use $f_s/2$.
The Mathworks docuementation for firpmord:
"Description
[n,fo,ao,w] = firpmord(f,a,dev) returns the approximate order n, normalized frequency band edges fo, frequency band amplitudes ao, and weights w that meet input specifications f, a, and dev.
example
[] = firpmord(,fs) specifies a sampling frequency fs. fs defaults to 2 Hz, implying a Nyquist frequency of 1 Hz. You can specify band edges scaled to a particular application's sample rate. You can use this with any of the previous input syntaxes."
In this case the maximum frequency that can be given in the f vector is 1.0 which corresponds to $f_s$. This is likely because the Matlab code is designing real (not complex) filters and you can only specify the response from 0 -> $f_s/2$ since the other half will be determined by conjugate symmetry.
It can be confusing - looking at Matlab's Remez and Remezord functions, they tend to say they normalize by the Nyquist frequency (they mean $f_s/2$) but they don't explicitly tell you that. Looking at the graphs of the filter responses in their examples in the documentation does make it clear - they show only half the frequency range and the maximum frequency is 1.0 (meaning it has been normalized by $f_s/2$.
In other code, e.g. ones that deal with complex filters, you may see the normalization use $f_s$, because you no longer have the conjugate symmetry you can specify the response over the whole frequency range (0-> $f_s$).
So it really depends. There are really two conventions for normalized frequency and you have to read the documentation or ask what whether they've normalized by $f_s$ or $f_s/2$.


coding gain and shaping gain in SCMA It is said that sparse code multiple access (SCMA) has the advantage of higher shaping and coding gain over LDS-CDMA due to the multidimensional constellation. What exactly is the shaping gain and how multidimensional constellation  increases the shaping gain in SCMA?

You can read about shaping gain in any book about digital communication. For example, John G. Proakis, "Digital Communications".
The theory is presented in J. H. Conway and N. J. A. Sloane, "Sphere Packings, Lattices and Groups".

In SCMA case, the multidimensional mother constellation only potential can increase shaping gain. There are many different methods for design good SCMA codebooks. I conducted a study comparing known codebooks and presented method for design codebooks based on genetic algorithm. For AWGN channel, there are codebooks based on QPSK alphabet with very good performance (in comparison with codebooks based on good multidimensional constellation).

In addition to a good multidimensional mother constellation, it is necessary to optimize the rotation factors, power variation coefficients, permutation, etc... Only a good multidimensional mother constellation (it has a shaping gain) does not guarantee the good performance of SCMA. So, QPSK mother constellation without shaping with good optimization of the rotation factors may provide the good performance of SCMA system.


"Tail-like" response from digital microphone I was doing some experiments with audio and wish to get "raw and unprocessed" data, and noticed that all the files had the following profile (just blowing into the mics):

Speaking to the microphone vendor, they mentioned - audio acoustic system will convert impulses into tail recovery dependent on the delay constants of the system which is dependent on system bandwidth.

Can someone explain what this means? What is the phenomenon called? Equations would be helpful as well.

Thanks

The tail is most likely the impulse response of a high pass filter. All acoustic system  are highpass by nature since air really can't transmit DC pressure. In addition many data acquisition systems have DC blocking filters and/or high passes build in to keep undesirable "rumble" signals that's below the range of hearing out of the recording.

Could also be in the microphone itself. It looks like you are hitting really hard and the recording is clipped. It takes a bit of time for the microphone and the pre-amp to recover from a major overdrive.


ANSI S1.11.2004 Octave-Band Filter implementation I am trying to implement some code to firstly generate an Octave-Band and Fractional-Octave Band filter.  I have used MATLAB's octaveFilter method to generate compliant Class-1 and Class-0, 1/1 and 1/3 Octave band filters, but I want to understand how to design these filters from scratch so I can better understand their construction.  Then I want to go on to implement in code.  So the coding environment doesn't matter, but I want to understand the filter construction from the filter definition/algorithm.
My problem is that the documentation for MATLAB's octaveFilter method doesn't explain the algorithm used to create the filters, nor does it explain how the filter is applied to data programmatically.
Can someone please point me to some material that explains how to create these filters, for any centre frequency I choose, and then implement in code?  Many thanks!

That's actually rather simple.
A little bit of trial and error will do the trick. The standard specifies upper and lower limit for the transfer function of the filter. So you can just try a normal Butterworth bandpass and see what happens.
The picture below shows a 3rd octave design at 1kHz and a sample rate of 48 kHz. A 2nd order bandpass isn't quite a good enough for class 1 limits, but a 3rd order does the job just fine. The 4th (and higher) order works as well.
Once you have a design that works at 1kHz, you can apply the same design to every other center frequency as well. At very high frequencies you will see some bilinear frequency distortion. Depending on your application you want to adjust for this or not.
Things to tweak are filter order, filter type and exact location of the corner frequencies. In the example the "nominal" corner frequencies worked just fine.


How to make semilogx to get a plot gain against log frequency with several values? [beginner] Gain: 20 log (Vout/Vin) with Vin = 1.0. I've been searching for gain values: a. -156.48 dB b. -80 dB c. -53.98 dB d. -13.98 dB e. 3.52 dB f. 10.32 dB g. 13 dB h. 13.93 dB i. 13.96 dB j. 13.98 dB

I am still at the very start on this semester to learn about matlab. How to make semilogx to get a plot gain against log frequency with several values? I've been searching matlab tutorial but I still got no clues. If anyone could give some explainations or examples relate with this question, please tell me. I want to learn about it. Thank you.


Is there a rule of thumb for zero padding in image processing? I see there are a lot of answers on why zero padding is necessary and how it avoids wrapping around the sides of images. However is there a rule of thumb on how much padding will be good for the image processing? I am mainly thinking from the point of view of:

a.) speed of processing the image

b.) noise

I have till now referred to the following answers:

Why images need to be padded before filtering in frequency domain
Advantages/disadvantage of zero padding
FFT Zero Padding - Amplitude Change

Response will be really appreciated.

Thank you.

See Press et al. "Numerical Recipes in C++". Chapter 13 on "Fourier and Spectral Applications", Section 13.1 with a subsection entitled "Treatment of end effects by zero padding".  This is perhaps the best summary of zero-padding anywhere in the literature (but it's not referenced in either the contents or index).  My edition of this book is the 2002 second edition and the page number there is 545.


Generation of a stream of data bits 0s and 1s at a certain bit or data rate This is a very basic question ,but basics important, : How is a stream of data bits (zeros and ones that express a certain information) and at a certain rate being generated? what is involved in terms of processing and possibly hardware?  a clock? a quantizer? a sampler? or what exactly?  I really would like to get a brief idea

Multimedia signals such as video,picture, music, speech have bandwidths associated with their baseband analog source nature.

For digital processing this bandlimited analog source is passed through a analog to digital converter (into a sequence of samples), which operates at a clock-frequency determined by the signal sampling rate, $F_s$, which is also limited by the source bandwidth due to Nyquist-Shannon sampling theorem.

Furthermore for each sample, a number of bits are used to quantize (numerically resolve) it. The number of bits is dictated by the signal amplitude accuracy requirement. For digital video, it's mostly 8-bits per RGB channel. Or for audio it can be 16-bits per stereo L/R channel.

Hence, the multiplication of samlping rate $F_s$ and bit-per-sample gives you the raw source bit rate $R_b$. You may prefer to apply some statistical or perceptual data compression to reduce $R_b$ to a minimum.

Then you choose a digital modulation scheme such as FSK,PSK,GMSK,OFDM... and transmit those bits through your favorite medium. Before doing this, you may increse $R_b$ (by adding redundancies) for channel transmission error correction purposes. The result is your bit-stream (not necessarily logic signal levels of 0s and 1s though) that's traveling accross the universe...

I hope I clarified it.


Verification of all steps needed to extract MFCCs from raw audio data I need some verification on steps required to extract MFCCs from raw audio as I'm finding that some sources outline them slightly differently.

Get raw audio frame
Run window function (eg Blackman) on raw audio data by multiplying each datapoint with the window
Perform FFT on windowed audio frame.
Calculate magnitudes by powering and summing real and imaginary parts of first half of FFT real^2 + imaginary^2
Convert magnitudes to dBv using 10 * Math.log10(magnitude)
Construct Mel Filter Bank consisting of eg 26 triangular filters spaced out across entire spectrum
Multiply dBv data with Mel Filter Bank to get 26 coefficients
Run log10 on each of the 26 coefficients
Perform DCT/FFT or ?iFFT? on 26 coefficients
Result produced is 26 MFCCs?

Regarding point 7: I'm not certain of input for multiplying with Mel Filter Bank. Should it be dBv, magnitudes, magnitudes converted to Mel scale or something else?

Point 9: I found some sources saying that said 26 coefficients should then be used as an input for DCT/FFTs as though they would be a signal. Other sources say the should be used for an inverse FFT as the result is mean to be a cepstrum. Can someone clarify?

Can someone verify this process and point out any mistakes?

I find it difficult to find solid points of reference covering these techniques and providing reasoning behind them. I was wondering if there's a book covering all of the above with technical examples, starting with theory and explanation of DSP concepts along with MFCCs and related speech analysis techniques. Most books I'm seeing focus mainly on Mixed Markov Models and don't seem cover cepstrums and MFCCs at all (at least by looking at their tables of contents).

Point 7: Multiply and add the FFT with the MelFilter Banks. Therefore, if you are using 26 Filter banks , we will end up getting 26 dimensional vector at the end.
Point 9: Take the DCT of the 26 dimensional vector.
Point 10: Keep the first 12 coefficients and discard the rest.

I understood MFCC features by reading this link: http://practicalcryptography.com/miscellaneous/machine-learning/guide-mel-frequency-cepstral-coefficients-mfccs/

It is very simple to understand from the link above.