Is initializing a digital filter's output with no "momentum" a non-trivial task? I'm working on implementing a filter with a very slow step response. This filter is implemented as a cascaded second-order-section filter (transposed direct form 2). I'm using the output of this filter as the input to a controller. Thus I'm trying to slow down how quickly the controller set point is able to change.

I'm running into an issue where I would like to initialize the set point of this controller. In order to do this, I need to provide initial conditions to the filter.

Let's assume I have a 7th order filter implemented in SOS stages as described above. Is it is a trivial task to find initial conditions such that I can specify the output but also have it such that the there is no "momentum" in the filter for it to drive from its initial output?

If this were an analog filter, I am essentially trying to specify that y(0) = constant, while y'(0), y"(0), ... etc. all equal zero.

I've never attempted to do something similar with a digital filter. Is this trivial? Is this non-trivial? Are there any references on the subject?

Thanks for the help!

Edit: So I've changed up my system design such my filter output isn't the setpoint, but rather a deviation from a setpoint. This allows me to simply zero the filter to achieve what I was originally trying to do. However, my question still stands for a matter of interest.

If you want a stable linear time-invariant system to output constant $y$, it must have received input $x$ that is the ratio of the constant output and the zero frequency response of the system $H(1):$

$$y = H(1)x\quad \Leftrightarrow \quad x = \frac{y}{H(1)}$$

Or, you'd want to change the state of the system so that it reflects that situation. The system must satisfy the condition $H(1) \ne 0.$

For a composite filter that consists of a serially connected cascade of second-order sections, you can start with the first section with input $x$ and calculate sequentially using $y = H(1)x$ (recycling the variable names to mean the inputs and outputs of a single section) the constant intermediate output of each section. Or you could start from the other end and use $x = y/H(1).$ Knowing both the input $x$ and output $y$, use the signal flow diagram of each section to calculate its state variables, starting with the summation point with dependency only to known constants. Because all oscillation has settled due to constant input and stability of the filter, all signals are constants rather than functions of time, and the delays in the flow diagram have become identities.

Figure 1. Transposed direct form II biquad signal flow diagram. (CC BY-SA 3.0 by Fcorthay)

Given constant input $x$, the constant output $y$ of a biquad section decipted in Fig. 1 is:

$$y = H(1)x = \frac{b_0+b_11^{-1}+b_21^{-2}}{1+a_11^{-1}+a_21^{-2}}x = \frac{b_0 + b_1 + b_2}{1 + a_1 + a_2}x$$

As for the state variables, the top sum $m_0,$ the middle sum $m_1,$ and the bottom sum $m_2$ are calculated as (start calculation from the bottom):

$$\begin{align}
m_0 &= b_0x + m_1\\
m_1 &= b_1x - a_1y + m_2\\
m_2 &= b_2x - a_2y\\
\end{align}$$

The top sum $m_0$ should equal $y$, which we can verify:

$$\begin{align}m_0 &= b_0x + b_1x - a_1y + m_2\\
&= b_0x + b_1x - a_1y + b_2x - a_2y\\
&= b_0x + b_1x - a_1\frac{b_0 + b_1 + b_2}{1 + a_1 + a_2}x + b_2x - a_2\frac{b_0 + b_1 + b_2}{1 + a_1 + a_2}x\\
&= \frac{b_0 + b_1 + b_2}{1 + a_1 + a_2}x\\
&= y
\end{align}$$


Matched filter - do I need to know the signal I am trying to find? I need to identify a matched filter, and so have gone through the mathematics to do this as listed here https://en.wikipedia.org/wiki/Matched_filter which derives the optimal filter $h$ to apply to a time series $x$, where that time series is assumed to comprise a desirable signal $s$ and some noise $v$. This suggests the optimal matched filter is:
\begin{equation}
    h = \frac{1}{\sqrt{s^{\mathrm{H}}\mathbf{\mathrm{R}}_v^{-1}s}} \mathbf{\mathrm{R}}_v^{-1}s
\end{equation}
and this should be applied to produce the filtered series $y$ thus:
\begin{equation}
    y\left(n\right) = \sum_{k=-\infty}^{\infty} h\left(n-k\right)x\left(k\right)
\end{equation}
however, I don't know $s$.  How, therefore, do I use this method (if at all)?

As stated by others, the signal $s[n]$ has to be known to calculate the corresponding matched filter.
An alternative is to use the power of deep learning to emulate a matched filter response as shown in this paper.


Power spectral density of $\left(x(t)\right)^2$? The relation between  $x(t)$ and output $y(t)$ of a non-linear device is expressed as
  $$y(t) = (x(t))^2$$
  Let $x(t)$ be zero-mean stationary Gaussian random process with auto-correlation
  $$R_x(\tau) = e^{-a\lvert \tau \rvert} , \quad a>0$$
  Find the output PSD?

I tried using the formula of output $\textrm{PSD}_y = \lvert H(f)\rvert^2 \cdot \textrm{PSD}_x$ , but stuck in process of finding $H(f)$ of a non linear system squarer. Not able to think of any other idea.

Since the question has been raised as to whether the hint that I had given to the OP in a comment on the original question was appropriate for a newcomer to signal processing, here goes.

Stripped of extraneous baggage and notation, the question is whether it is possible to determine the value of $E[X^2Y^2]$ straightforwardly where   $X$ and $Y$ are zero-mean unit variance jointly Gaussian random variables with correlation coefficient $\rho$, or is it necessary to resort
to advanced methods or results such as the one cited in the answer by @MattL or generalizations thereof such as Price's Theorem (cf. the
book cited in MattL's answer).  One straightforward method involves
conditional expectations and the result that is known to be a LIE:

$$E[X^2Y^2] = E\big[E[X^2Y^2\mid Y]\big]$$

Sorry, I couldn't resist: LIE is an acronym for the Law of Iterated Expectation which asserts that the expectation of $A$ is the
same as the expectation of the random variable $E[A\mid B]$. Note that $E[A\mid B]$ is a function of $B$, not $A$ as one might naively believe. We already squeezed out all of the randomness in $A$ when we took its
expectation. The LIE is
that the expected value of $E[A\mid B]$, a function of $B$) works out, by a miracle of
modern mathematics, to be the same as the expected value of $A$.

Now, for any choice of real number $y$, given that $Y=y$,

$$E[X^2Y^2\mid Y = y] = E[X^2y^2\mid Y = y]
= y^2 E[X^2\mid Y = y].$$
But, given that $Y = y$, we know that the conditional distribution
of $X$ given $Y=y$ is a Gaussian distribution with mean $\rho y$
and variance $1-\rho^2$. Consequently,
$$E[X^2\mid Y = y] = \operatorname{var}(X\mid Y=y)
+ \left(E[X\mid Y = y]\right)^2 = (1-\rho^2) + (\rho y)^2$$
leading to
$$E[X^2Y^2\mid Y = y] = y^2((1-\rho^2) + (\rho y)^2))
= y^2(1-\rho^2) + \rho^2y^4 $$
giving us that the random variable $E[X^2Y^2\mid Y]$ is just
$Y^2(1-\rho^2) + \rho^2Y^4$. Since $E[Y^2] = 1$ and $E[Y^4] =  3$,
we get that
$$E[X^2Y^2] = E\big[E[X^2Y^2\mid Y]\big] = E[Y^2(1-\rho^2) + \rho^2Y^4]
= (1-\rho^2) + 3\rho^2 = 1 +2\rho^2.$$

Finally, noting that $\rho = e^{-a|t|}$, we can write
$$R_y(t) = 1 + 2e^{-2a|t|}$$
just as has been explained in MattL's answer.


How to choose the number of codewords in bag of words model In the bag of words model, generally is it better to have more codewords or less? For example, if I have 1000 x 20 raw features, then what would be good value for codebooks? I am thinking of smaller size codebook ( such as 20) so that it forces the representation to be sparse.

It is a trade-off between recall and precision. If you have only one codeword in your codebook then everything will map to it and you will have high recall. If you have 20,000 codewords in your codebook then you will have high precision and low recall. You can use a validation set to try different codebook sizes and see which one works best for your application.

If you have a smaller size codebook I don't think that will help you get a sparse feature. It will just mean that more features map onto the same small set of codewords, increasing the recall.


Derivation of the LMMSE (Linear Minimum Mean Squared Error) Estimate and the MMSE Under Gaussian Prior I am learning estimation theory through Steven M. Kay - Fundamentals of Statistical Signal Processing, Volume 1: Estimation Theory. In the Chapter 12 (Linear Bayesian Estimator), Theorem 12.1 (Bayesian Gauss-Markov Theorem) gives the LMMSE estimation of the signal based on the linear noisy measurement under the Gaussian prior assumption:

If the data are described by the Bayesian linear model form
\begin{equation}
\boldsymbol{x}=\boldsymbol{H\theta}+\boldsymbol{w} \tag{12.25}
\end{equation}
where $\boldsymbol{x}$ is an $N \times 1$ data vector, $\boldsymbol{H}$ is a known $N\times p$ observation matrix, $\boldsymbol{\theta}$ is a $p \times 1$ random vector of parameters whose realization is to be estimated and has mean $E(\boldsymbol{\theta})$ and covariance matrix $\boldsymbol{C}_{\theta\theta}$, and $\boldsymbol{w}$ is an $N \times 1$ random vector with zero mean and covariance matrix $\boldsymbol{C}_w$ and is uncorrelated with $\boldsymbol{\theta}$ (the joint PDF $p(\boldsymbol{w},\boldsymbol{\theta})$ is otherwise arbitrary), then the LMMSE estimator of $\boldsymbol{\theta}$ is
\begin{align}
\hat{\boldsymbol{\theta}} & = E(\boldsymbol{\theta})+\boldsymbol{C}_{\theta\theta}\boldsymbol{H}^T(\boldsymbol{H}\boldsymbol{C}_{\theta\theta}\boldsymbol{H}^T+\boldsymbol{C}_w)^{-1}(\boldsymbol{x}-\boldsymbol{H}E(\boldsymbol{\theta})) \tag{12.26} \\
& = E(\boldsymbol{\theta})+(\boldsymbol{C}_{\theta\theta}^{-1}+\boldsymbol{H}^T\boldsymbol{C}_w^{-1}\boldsymbol{H})^{-1}\boldsymbol{H}^T\boldsymbol{C}_w^{-1}(\boldsymbol{x}-\boldsymbol{H}E(\boldsymbol{\theta})) \tag{12.27}
\end{align}
The performance of the estimatior is measured by the error $\boldsymbol{\epsilon}=\boldsymbol{\theta}-\hat{\boldsymbol{\theta}}$ whose mean is zero and whose covariance matrix is
\begin{align}
\boldsymbol{C}_\boldsymbol{\epsilon} &= E_{\boldsymbol{x},\boldsymbol{\theta}}(\boldsymbol{\epsilon}\boldsymbol{\epsilon}^T) \\
& = \boldsymbol{C}_{\theta\theta} - \boldsymbol{C}_{\theta\theta}\boldsymbol{H}^T(\boldsymbol{H}\boldsymbol{C}_{\theta\theta}\boldsymbol{H}^T+\boldsymbol{C}_w)^{-1}\boldsymbol{H}\boldsymbol{C}_{\theta\theta} \tag{12.28} \\
& = (\boldsymbol{C}_{\theta\theta}^{-1}+\boldsymbol{H}^T\boldsymbol{C}_w^{-1}\boldsymbol{H})^{-1} \tag{12.29}
\end{align}

Since the prior of $\boldsymbol{\theta}$ is Gaussian, the LMMSE estimate $\hat{\boldsymbol{\theta}}_{LMMSE}$ is equivalent to the MMSE estimate $\hat{\boldsymbol{\theta}}_{MMSE}$, and $\hat{\boldsymbol{\theta}}_{MMSE}$ is equal to the posterior mearn $E(\boldsymbol{\theta}|\boldsymbol{x})$. Since the prior and likelihood are both Gaussian, the posterior distribution $p(\boldsymbol{\theta}|\boldsymbol{x})$ is also Gaussian.
Here I am trying to derive $\hat{\boldsymbol{\theta}}_{MMSE}$ and $\boldsymbol{C}_\boldsymbol{\epsilon}$ from the perspective of PDF multiplication, that is, calculate $p(\boldsymbol{\theta}|\boldsymbol{x}) \propto p(\boldsymbol{x}|\boldsymbol{\theta})p(\boldsymbol{\theta})=\mathcal{N}(\boldsymbol{x};\boldsymbol{H\theta},\boldsymbol{C}_{w})\mathcal{N}(\boldsymbol{\theta};E(\boldsymbol{\theta}),\boldsymbol{C}_{\theta\theta})$, and formulate the quadratic and firse-order terms of $\boldsymbol{\theta}$ at the exponential to form a Gaussian PDF. The covariance matrix of $p(\boldsymbol{\theta}|\boldsymbol{x})$ I got matches 12.29, but the posterior mean is the following form:
\begin{equation}
E(\boldsymbol{\theta}|\boldsymbol{x}) = \boldsymbol{C}_{\boldsymbol{\epsilon}}(\boldsymbol{H}^T\boldsymbol{C}_w^{-1} \boldsymbol{x}+\boldsymbol{C}_{\theta\theta}^{-1}E(\boldsymbol{\theta})) \tag{q1}
\end{equation}
So my question is, is the posterior mean I got in q1 equal to the $\hat{\boldsymbol{\theta}}$ given in 12.26 and 12.27? If so, how can I reach that?
By the way, I can't find the way from 12.26 to 12.27 (12.28 to 12.29 either). So can someone give me a hint?

In the past I derived it as following:

It is a little different approach.
If it answers your question I will rewrite it in a proper LaTeX.
Regarding your question about the steps in the derivation you presented, it is using the Woodbury Matrix Identity (Both 12.26 to 12.27 and 12.28 to 12.29).
Related Answers:

Deriving the Matrix Inversion Lemma for RLS Equations vs the Woodbury Derivation.
Sequential Form of the Least Squares Estimator for Linear Least Squares Model.

In the answer above you may see some related derivations.


How to remove constant accelerometer velocity offsetting? For quite some time I have been working on a project and got stuck recently with a problem considering integrating vertical acceleration from a legendary MPU6050 to get vertical velocity.
What I have been struggling with is the fact that velocity does not ever come back to 0. After every sequence of movement it just gets more and more off. It doesn't drift linearly anywhere, just stays constant, but wrong. It kind of seems, like the accelerometer was more sensitive detecting the acceleration in one direction than in the opposite one.
What may cause the problem? Is there any way to fix it with software? Switching to an IMU like BNO055, that is capable of doing the integration maths 'in chip' isn't a possibility right now. Thanks for all the help in advance.
Have a nice day.

I guess you are sure the device is definitely at rest (v=0) after each event, and not just moving slowly?
If so a couple of possible causes to consider:
NUMERICAL: You might get this type of thing if the integration is inaccurate.
Consider reducing sample interval, and use higher order integration.  However, I would
have expected numerical drift to be sometimes positive, sometimes negative. If it is always negative, numerical seems less likely to me.
INSTRUMENTAL: Does that accelerometer need to be "level"? Is it possible to repeat
the experiment with the accelerometer mounted upside down to see
if the effect reverses?
An immediate suggestion to better understand the problem (which you have probably
already done) is compare the acceleration and velocity plots. As I recall once
the system is underway the two curves should be 90 degrees out of phase, so there
should be time alignment between peaks and zero crossing as in this sketch.

You could check the time alignment at those points. It might indicate if the drift
is arising at a particular part of the cycle (as you suggest).
Pragmatic Adjustment:
A perfect correction would require identification of the exact cause. However, depending on what you are using the velocity curve for it might be acceptable to implement a "pragmatic" fix in processing. As one idea, in an ideal world the areas under the two positive peaks (in accel plot) should balance the area above the negative peak.  The observed drift (dv in picture) should equate to this difference. Hence for example, you may be able to introduce a percentage boost
to the positive acceleration values. With a little thought this should be able
to be automated.
However first I'd suggest looking at the time alignment and see if that reveals anything useful.


Phase shift and phase spectrum terms in multidimensional signal I know about phase of a 1D signal. But when I go into higher dimensions like 2D,3D etc, it becomes headache to grasp the concept.

What are the terms phase shift and phase spectrum mean in case of multidimensional(2D,3D etc.) signal like an Image or video signal?
How is the phase term represented mathematically?
How do we analyze phase spectrum?
Which information do we get from phase spectrum about the signal?
Can anybody explain with a plot or figure ?

For a real signal, the magnitudes of the Fourier transform bins tell what are the amplitudes of the sinusoidal components that an infinite tiling of the signal frame can be decomposed into. The phases of the bins encode translation of the sinusoids along their direction, like so:

The translation can be decomposed into a series of component translations one along each of the axes. The axes of an image are spatial axes, and in video the third axis is temporal. The component translations have magnitudes inversely proportional to the frequencies indicated by the multidimensional index of the bin. For example, you can easily spatially translate an image in the frequency domain by adding to all phases (in a polar form representation of the bins) a function of form $af_x + bf_y$ where $a$ and $b$ are constants that give the amount of translation and $f_x$ and $f_y$ are the vertical and horizontal frequencies (bin indexes).

Translation of a sinusoid of zero frequency is a bit of an oxymoron. With real FFT its phase is fixed to zero.


Segmenting a sound file before each of its attacks Like the title says...
I want to split audio file where a single instrument is played (not so fast)...
I just need to detect the attack phase of the waveform to know where a played note begins. Can someone suggest a way to do this? It should be robust, not to react to isolated peaks...
If there are open source implementations please share.... or even pseudocode is fine ?
Cheers

If it is a single instrument playing without any background noise/instruments, then the problem becomes simple. All you have to do is keep track of energy, which is calculated at a defined window size (say 40 milliseconds). Wherever you see the energy increasing and sustaining for a while (you'll have to define how long a note is good for you to split your audio), you can mark the ascent of such an energy plot as the point to split.

You can also split audio when a particular note is played, for this you'll have to track energy in the respective spectrum band of the note you want to split the audio on. For example to track note A4 (a') which has 440Hz frequency. You might have to track the energy in a band around 440 Hz, say 420 to 460 Hz. Mind you, this splitting business based on notes played gets more complicated when complex notes are played.

Also, this problem gets really difficult when you have multiple instruments.


Fourier transform of a normalized vector I was expecting that when I take (for example, in MATLAB), s2 = fft(s1/norm(s1)), the norm of s2 is not equal one, norm(s2) $\neq$ 1.

I cannot figure out why this is! Can anyone shed light on this?

There are different conventions in scaling of the FFT, in MATLAB you need to scale it by $\sqrt{N}$, where $N$ is your number of samples. Saying it in matlabish:

clc, clear all
%% Create some sinusoidal signal with noise
t = linspace(0, 6*pi, 1000)+randn(1,1000);
s1 = sin(2*pi*t);
% Calculate the norm of s1
s1_norm = norm(s1);
display(sprintf('L2 norm of s1: %.2f', s1_norm))
%% Create signal s2 as the FFT of normalised s1
% Scale FFT by sqrt(N) - this is convention used in MATLAB
s2 = fft(s1)/sqrt(length(t));
% Calculate the norm of s2
s2_norm = norm(s2);
display(sprintf('L2 norm of s2: %.2f', s2_norm))


Estimating the Signal by Deconvolution with a Prior on the Filter Coefficients and the Signal Samples Assume I have signal $y[n]$ which is a result of convolution between channel $h[n]$ and signal $x[n]$. which means:

$$y[n] = h[n] \ast x[n]$$ where $\ast$ is the convolution operation

The signal $y[n]$ could be complex since we can consider the channel $h[n]$ to be complex too.

In normal case, $h[n]$ should be known in order to estimate $x[n]$ using deconvolution process.

My question what about if I only know that $ \forall n$, $$x[n] \in \left\{ -1, 1 \right\}$$ it means $x[n]$ is either $1$ or $-1$, it's a vector and each value in the vector is $1$ or $-1$ ($x[n]$ is a vector of +1/-1 we means I need to estimate the whole vector).

For example if its length is $4$ , it could be $[1,-1,-1,1]$ or $[1,1,1,1]$ and so on. Is it possible in that case to perform anyway, for example, deconvolution or any other method to estimate $x[n]$?

NP: The maximum length of $x[n]$ can be $256$ and maximum length of $h[n]$ can be $64$ , and we have a known information about vector $x[n]$ which is the $sum(x[n]) = 0$.

I would take approach based on Blind Deconvolution.

Since we're dealing with ill posed problem some assumptions should be made.
The intuitive approach would be using the information as a prior for the signal. Another idea is to add LPF assumption of the Filter by setting the sum of its coefficients to be 1 and non negative. Yet since we have Discrete Prior on the signal we're getting into combinatoric problem.
Which means brute force solution where the number of combinations is $ {2}^{n} $ where $ d $ is the number of the signal samples.

For $ n \leq 16 $ I'd say it would work for the given input size.
Yet for solution with higher number of samples this method isn't feasible.

In order to deal with higher dimensions (More samples) I'd use a GMM:

Namely the Prior model is a Gaussian Mixture Model (GMM) with 2 Gaussian's centered at $ \left\{ -1, 1 \right\} $ with very small variance in order to approximate discrete probability function.

So the problem I'm looking to solve is given by:

$$\begin{aligned}
\arg \min_{h, x} \quad & \frac{1}{2} {\left\| h \ast x - y \right\|}_{2}^{2} \\
\text{subject to} \quad & \sum h = 1, \, h \succeq 0
\end{aligned}$$

Step 1 - Solving for the Filter $ h $

Given the signal $ x $ is known, solving for the filter is pretty easy using matrix form of the problem (Which is convex):

$$\begin{aligned}
\arg \min_{h} \quad & \frac{1}{2} {\left\| X h - y \right\|}_{2}^{2} \\
\text{subject to} \quad & \sum h = 1, \, h \succeq 0
\end{aligned}$$

I showed, code included, how to solve such problem in my answer to How to Project onto the Unit Simplex as Intersection of Two Sets (Optimizing a Convex Function)?

Step 2 - Solving for the Signal $ x $

The model is $ y \mid h \sim \mathcal{N} \left( h \ast x, {\sigma}_{n} I \right) $ and the prior $ {x}_{i} \sim 0.5 \mathcal{N} \left( {\mu}_{1} = -1, {\sigma}_{1}^{2} = {0.1}^{2} \right) + 0.5 \mathcal{N} \left( {\mu}_{2} = 1, {\sigma}_{2}^{2} = {0.1}^{2} \right) $.

I'd use the MAP so we have:

$$\begin{aligned}
\arg \max_{x} p \left( x \mid y \right) & = \arg \max_{x} p \left( y \mid x \right) p \left( x \right) \\
& = \arg \max_{x} \log p \left( y \mid x \right) + \log p \left( x \right) \\
& = \arg \min_{x} -\log p \left( y \mid x \right) - \log p \left( x \right) \\
& = \arg \min_{x} \frac{1}{2} {\left\| h \ast x - y \right\|}_{2}^{2} - \lambda \log p \left( x \right) \\
& = \arg \min_{x} \frac{1}{2} {\left\| h \ast x - y \right\|}_{2}^{2} - \lambda \sum \log p \left( {x}_{i} \right)
\end{aligned}$$

Where $ \lambda \propto N {\sigma}_{n}^{2} $ where $ N $ is the number of samples (Dimension of $ y $).

This can be solved by any solver. I'd use MATLAB's fminunc().
Though one could alter (For MAP Estimation) the Expectation Maximization (EM) process for faster and better converging algorithm. Another option would be using Probabilistic Programming with one of the options available today.

Remark: The above is Bayesian Modeling of the problem. One could build optimization problem with some intuition in the form of:

$$ \arg \min_{h, x} \frac{1}{2} {\left\| h \ast x - y \right\|}_{2}^{2} + \lambda \sum_{i = 1}^{m} {\left( {x}_{i}^{2} - 1 \right)}^{2} $$

Which isn't motivated by Bayesian model but still drives the solution to where we want it.

Initialization

One approach to initialization of the estimated $ x $ signal is to use hard threshold. So we set $ {x}_{i} = 1 $ if $ {y}_{i} \geq 0 $ and $ {x}_{i} = -1 $ otherwise.

This approach could also be a greedy method to solve the step for $ x $. Yet it doesn't take under account the delay of the filter.

After iterating enough for a stable solution one could round the result of $ x $ such that $ {x}_{i} \in \left\{ -1, 1 \right\} $.

I haven't tested this approach myself, but I really like its model.
I will publish the results of MATLAB simulation soon.


Verifying Channel Response Identity If $x[n]$ is input signal, $y[m]$ is output signal, $h[m]$ is impulse response of the channel, then the following relation holds

$\lvert Y(\omega)\rvert = \lvert H(\omega)\rvert^2 . \lvert X(\omega)\rvert$

where $F(\omega)$ is the DFT of $f[n]$.

To verify this, I wrote following MATLAB code

x = randn(1,1000);
h = [1 2 3 4 5];
y = conv(x,h);
plot((abs(fft(h,1024)).^2).*(abs(fft(x,1024))));
hold on
plot(abs(fft(y,1024)),'r')

Result shows that second plot is at a much lower level than first plot. When I multiply second plot by $1024$ i.e. NFFT then both the plots overlap.

What is wrong with the original equation? I remember learning it that way in class.

I think you made a mistake between the power spectrum and the Fourier Transform.
This is the right form of the Convolution Theory (Nothing is squared):

$$ \lvert Y(\omega)\rvert = \lvert H(\omega) X(\omega)\rvert $$

Try this and it will work for you on MATLAB.

MATLAB Code

% Calculating the Magnitude of the Filtered Signal

vFilterCoeff = [1; 2; 3; 4; 5];
vInputSignal = randn(1000, 1);

vOutputSginal = conv(vInputSignal, vFilterCoeff);

fftNumPoints = 1024;

vOutputSignalDftMagnitude = abs(fft(vOutputSginal, fftNumPoints));
vFilteredSignalDftMagnitude = abs(fft(vFilterCoeff, fftNumPoints) .* fft(vInputSignal, fftNumPoints));

hFigure = figure();
hAxes   = axes();
plot([vOutputSignalDftMagnitude, vFilteredSignalDftMagnitude]);

The result is:

As can be seen there's a full equality.


FIR filter design: evaluating the error function Question 1

$$
H(e^{j\omega})=\sum_{n=0}^{N-1}h[n]e^{-jn\omega}
=\mathbf{c}^H(\omega)\cdot \mathbf{h} \tag{1}
$$
$$
=\mathbf{h}^H\cdot\mathbf{c}(\omega)  \tag{2}
$$

$$H(\mathbf{h})=\sum_{k=1}^Kh[k]e^{j\omega_k}\tag{3}$$

The design of filter in matlab ,if using:

$-$$(1)$ and $(3)$ the length of filter is $K$.

$-$$(2)$ the length of filter is $N$.

What is the difference between these two design?
Can i use the operation 'transpose' in $D(w)$ and weighting function,length of filter?

Question 2

application in matlab

I used the same functions in 'cfirls' like desired and actual function.

   %% Desired coefficient

   N = 61;              %desired filter length
   tau = 26;            % desired passband group delay
   % frequency grid
   f = [linspace(-1,-.18,164),linspace(-.1,.3,80),linspace(.38,1,124)];
   % desired frequency response
   d = [zeros(1,164),ones(1,80),zeros(1,124)].*exp(-1i*pi*f*tau);
   w = [10*ones(1,164),ones(1,80),10*ones(1,124)];
   A = w(:,ones(1,N)) .* exp(-1i*pi*f*(0:N-1));
   h = A \ (w.*d);
   %%actual function
   c=[zeros(1,164),ones(1,80),zeros(1,124)]*exp(-1i*pi*f*(0:N-1));
   H=h'.*c;

   %% error function
   E=abs(H-d)^2;
   S=sum(E');

   plot(f,S)

I need in my work to reduce the error as little as possible.How can i changing in this code?

There are many problems in your question and in your code. First, Eqs $(2)$ and $(3)$ are wrong. Eq. $(2)$ cannot equal Eq. $(1)$ because it is its complex conjugate (and $\mathbf{c}$ is a complex vector). Eq. $(3)$ is meaningless.

I will not rewrite your Matlab code, but you should have added that it doesn't run.
There are many errors, among which several incompatible matrix dimensions. You should use the function cfirls instead of just copying some of its code. Then vector dimensions would be checked and care would be taken that all vectors have the same format (I chose columns).

Assuming that all these problems are resolved (simply by calling cfirls with the vectors f, d, and w), you can compute the error, but this is where the next problem occurs. Why don't you just use Eq. $(1)$ in your question to compute the frequency response? (You might as well use the Matlab command freqz). The way you define the matrix c in order to compute the frequency response H doesn't make sense.

Assuming that you manage to compute the designed frequency response on a grid, you could choose to compute the unweighted error, as suggested by your code (even though you minimized a weighted error). However, your variable S will be a scalar (as soon as the code to compute it is fixed). So how do you want to plot it as a function of frequency (as suggested by the last line of your code)?

In sum, apart from using cfirls, which will give you the correct filter coefficients, think about how you should compute the frequency response, and why you defined the matrix c the way you did. Then think about which quantity you would like to plot. As mentioned before, the squared error is just a number, not a function of frequency. It could be useful to plot the magnitude of the difference between desired and designed frequency response:

$$\big|E(\omega)\big|=\big|H(\omega)-D(\omega)\big|\tag{1}$$

With N, f, d, and w defined as in your code, the most straightforward way to do what you're trying to do would be:

h = cfirls(N,f,d,w);
H = freqz(h,1,pi*f);
E = H - d;
plot(f,20*log10(abs(E)))


Simple FFT convolution with impulse response wav file doesn't give the right results I have been trying to find a solution to my FFT convolution problem but everything that I have found while searching hasn't been helping me.

I have a working forward and backward FFT function, that I have tested with perfect results. A wav file that I ran through the forward FFT came out with a result of a frequency domain data, ran this domain data backwards again and back came a nice sounding waveform.

Here's the problem; I have now tried to put forward FFT to a pre-recorded piano sound (wav) a file with only one note and no dead silence. Then I've done the same to an IR-file (wav) of a room reverb, took the results from these two forward FFTs and tried to convolve them with the following code (the res_ir buffer has been scaled to have values between -1 to +1 before this):

long bufSize = 1024*62;
mFFT->runFFT(res, bufSize);
mFFT->runFFT(res_ir, bufSize);

double a=0;
double b=0;
double c=0;
double d=0;
for (int i=0;i<bufSize;i++)
{
     a=res[i].re();
     b=res[i].im();
     c=res_ir[i].re();
     d=res_ir[i].im();

    res[i].mRe = (a*c - b*d) * (gain);
    res[i].mIm = (a*d + b*c) * (gain);
}
mFFT->runIFFT(res, bufSize,true);

I'm not looping through the complete wave file but just the size of one really big FFT buffer (1024*62), this to make it easier for me to follow and to eliminate extra possible bugs or errors keeping the problem down in size.

the result I'm getting is just the IR wav sound shaped somewhat like the piano file, but it sounds like the "starter gun" from the IR file (no piano sound at all).

Here's an image if the results (starting from above, first is the IR wav with the starter gun, second is the piano wav and last is the result from this):

Someone who can help me find what or where to dig for answer to this problem?

Your code performs the circular convolution.
Circ Convolution gives an array equal to max( len(a) , len(b) )
Convolution will return an array of length = len(a) + len(b) - 1 // a, b > 0
If you think about it, adding a reverb makes a sound longer. That's the equation ^^ that describes by how much. What you're hearing is the end of the reverb played over the start.
Keep in mind in a musical setting, often the original sound is mixed in with the reverb sound.

I zero padding your data before FFTing it.

long bufSize = 1024*62;
for(int i=bufSize; i<2*bufSize; ++i) res[i].mRe = res[i].mIm = res_ir[i].mRe = res_ir[i].mIm = 0;
mFFT->runFFT(res, bufSize);
mFFT->runFFT(res_ir, bufSize);

double a=0;
double b=0;
double c=0;
double d=0;
for (int i=0;i<bufSize*2;i++)
{
    a=res[i].re();
    b=res[i].im();
    c=res_ir[i].re();
    d=res_ir[i].im();

    res[i].mRe = (a*c - b*d) * (gain);
    res[i].mIm = 0;
}

That gets rid of the 'wrapping' effect of circular convolution at least.
Also, the imaginary part is zero because the convolution of two real signals is real.


How can i define the vector to variable scalar? I have A,B,C,D,E values all in size 1x1155.
A=[1.0001, 1.0005, 1.00010,...,1.0009]. similar to B, C, D and E.
I need (A(1)+B(1))/(C(1)+D(1)E(1)) + (A(2)+B(2))/(C(2)+D(2)E(2)) + ... + (A(1155)+B(1155))/(C(1155)+D(1155)E(1155)). Therefore i need each individual variable using loop rather than defining one by one from 1 to 1155. Is there a way to do that in Matlab?

Use pointwise vector operations:

% given A,B,C,D,E
A = rand(1,1155);
B = rand(1,1155);
C = rand(1,1155);
D = rand(1,1155);
E = rand(1,1155);
% compute F
F = (A+B)./(C + D.*E);


For a system to be causal, number of finite zeros <= number of finite poles. Why? I read in this pdf that for a system to be causal, the number of finite zeros must be no greater than number of finite poles. Why?

I know that for a system to be causal, $h[n]=0$ for all $n<0$. Can you give an explanation how the statement maps to this? Or why the number of finite zeros are less than or equal to the number of finite poles for a system to be causal?

If the number of finite zeros is not greater than the number of finite poles then the transfer function is proper, i.e., the degree of the numerator polynomial is not greater than the degree of the denominator polynomial.

If the degree of the numerator polynomial were greater than the degree of the denominator polynomial, we would get at least one pole at infinity. This violates the requirement that a causal and stable system must have all its poles in the left half-plane (for continuous-time systems), or inside the unit circle (for discrete-time systems). Consequently, a causal and stable system must have a proper transfer function. This is true in continuous time as well as in discrete time.

These two answers are related: answer 1, answer 2


What Are the Types of Deconvolution? I am totally new to image processing and wanted to ask you if you could confirm what I understood. It is about deconvolution:

From what I read we find 2 main types of deconvolution:

1. Analytical

   i.Linear : inverse filter, Wiener deconvolution..
   ii. non linear : maximum  entropy  method , CLEAN algorithm..

2. numerical

Is this right ? If I need deconvolution for deblurring, how do I proceed to choose the right algorithm?
Thank you in advance!

I would say you can classify using the following main properties:

Blind Deconvolution.
Non Blind Deconvolution.

Then I'd follow:

Linear Model.

Time / Spatial Invariant Model.
Time / Spatial Variant Model.

Non Linear Model.

Time / Spatial Invariant Model.
Time / Spatial Variant Model.

I think those are the main properties of a Deconvolution Problem.
The hardest would be Blind Deconvolution with Non Linear Model of Variant Operator.


Analytic solution for non-flat filter design I'm trying to write an accelerometer calibration script that uses filters to convert from volts into $m/s^2$.  As accelerometers tend to have non-flat response curves, this means I have to design a rather complex filter.  I'm not worried about phase, as I can just apply the filter twice in opposing directions to correct for any phase offsets (like matlab's filtfilt), so the focus is on designing a filter that approximates a user-provided magnitude curve.

Ideally, the user provides a calibration curve as input into an analytic algorithm to solve for the best fitting filter poles.

I'm aware MATLAB has a filter design function, but I don't know what the underlying algorithm is (if its an optimizer, or a closed form solution).

So my question is...

Is there an analytic solution to my filter design problem?  Or do I have to use optimisation scripts to get the best filter?

I'm not mentioning programming language here, as I want to understand the underlying math behind this.

If you are not concerned on the phase and just want to approximate a magnitude response, then your first option should be the frequency sampling method implemented in Matlab/Octave fir2() function.

You would provide the frequency grid and corresponding frequency response magnitude at those frequencies.

As you have also mentioned, least-squares approach is another alternative. Indeed by using suitable weights, you can distribute the error according to your priority cirteria.

Magnitude approximation based on LMS adaptive system identification is also a possible option.


Are Computer Stereo Vision techniques appropriate for sub-millimetre measurements? I have a project where I would like to image an object and be able to derive the heights of features in this image to sub-millimetre precision (exactly how precise is still yet to be determined, but let's say 100ths of a millimetre for now).
I have been previously advised that direct laser ranging techniques will not be appropriate

the travel time will be too small and thus will require too much precision to make precise calculations
minor vibrations (such as a person walking near the apparatus) will perturb the results

I have observed a laser device that sells for approximately $1000 that can achieve the precision but suffers from the vibration problem (which is fine, mechanically isolating the apparatus is another discussion).
I would prefer to achieve a result that is more cost effective, and considered stereo vision as an alternative.  Being a novice to this field I am uncertain if the desired precision can be achieved.
Is the desired precision (at least) theoretically attainable?
Is there a recommended paper or resource that would help explain this topic further?
Additional notes
The objects in question will range from approximately 1/2" square up to about 2 1/2" square with some times very low thickness (1/16"?). A large percentage of the surface should be flat, though one test will be to confirm that assertion. Features will be fairly rough (generally sharp transitions). Aug 17 at 11:00
One of the "harder" interesting objects would be about 20mm square, 1.25mm high. The surface features in question would be on the order of .1 - .3mm I'm estimating. The camera position would likely be on the order of 6" above. Does this give you better insight? Aug 17 at 15:15
I am not looking to perform a single profile/relief measurement, but rather attempting to generate a surface height map of the object.  The surface features of the object, as well as the overall profile, are of significant interest.

Stereo imaging

Given the large field of view you need in relation to the accuracy you want, and how close you want to be, I think that stereo imaging may be a challenging, so you need to somehow amplify the differences you are trying to measure.

Structured lighting

If you are essentially trying to measure the profile of an object, have you considered a single high resolution camera and structured lighting?

Thanks to looptechnology for this image, used without permission, but hopefully attribution will be enough.

Note, the shallower the grazing angle, the the greater accuracy you can measure, but the lower the supported depth of field would be, so for your application you would need to optimise for your needs or make your system adjustable (one laser angle for 0-500um, another for 500-1500um and so on). In this case though, you would probably have to calibrate each time you changed laser position.

Incidentally, a very cheap way to try this out would be to pick up a pair of Laser Scissors which include a basic line laser LED.

Finally, you can remove the vibration problem by sampling multiple times, rejecting outliers and then averaging. A better solution though would be to mount the whole test apparatus on a block of granite. This worked well for laser micro-machining tools I've worked with in the past, which require micron level position and depth of focus accuracy, even when located in factories.

Some back of the envelope calculations.

Lets assume an incident angle of 10 degree from horizontal, and a camera with a 640x480 resolution and a field of view of 87 x 65mm. If we place the beam so that it is right at the bottom of the portrait frame with no sample, and then place the sample with the beam crossing it, this should give us a maximum height of around 15mm and thus an uncorrected resolution of around 24um for each pixel the line walks up the screen. With this setup, a 0.1mm variation should be visible as a 4 pixel variation in position.

Similarly, if we use an incident angle of 2 degrees from horizontal then this should give us a maximum height of around 3mm (Tan(2deg)*87mm) and thus an uncorrected resolution of around 4.7um per pixel, for a much more noticeable 20 pixel jump. This would probably require a much more accurate line laser however.

Note, if the camera is close enough then the you may need to do a second trig calculation, using the camera height, to determine the the true position of the line relative to the base line.

Also note that if you don't need absolute accuracy, and local repeatability is enough (say you are profiling the flatness of a sample to ensure it is within given tolerances) then just being able to see the relative position of the laser line might be enough.


Attenuate the low frequencies in a speech wave I am designing an FIR filter to attenuate the frequencies below 300Hz. I designed the filter using the scipy functions and the code is as below.

from numpy import cos, sin, pi, absolute, arange
from scipy.signal import kaiserord, lfilter, firwin, freqz
from pylab import figure, clf, plot, xlabel, ylabel, xlim, ylim, title, grid,
axes, show

nyq_rate = 48000 / 2
width = 100 / nyq_rate
ripple_db = 60.0
cutoff_hz = 300.0

N, beta = kaiserord(ripple_db, width)
taps = firwin(N, cutoff_hz/nyq_rate, window=('kaiser', beta), pass_zero='highpass')
w, h = freqz(taps, worN=512)

plot((w/pi)*nyq_rate, absolute(h), linewidth=2)
xlabel('Frequency (Hz)')
ylabel('Gain')
title('Frequency Response')
ylim(-0.05, 1.05)
grid(True)

show()

When I a execute it, the following error pops up.

ValueError: A filter with an even number of coefficients must have zero response at the Nyquist frequency.

I am not sure what is going wrong. Any help is appreciated.

As you've already figured out, you need to choose an odd filter length. I'll explain why. You're designing a linear phase FIR filter, and such filters satisfy certain symmetry conditions. See this answer explaining the four types of linear phase FIR filters (even/odd length, even/odd symmetry).

For frequency selective filters, such as the one you're trying to design, you need filters with even symmetry (type I or type II). An even length (type II) filter necessarily has a zero at Nyquist, so you can't use it to design a high pass filter. A type I filter (even symmetry, odd filter length) doesn't have this restriction, so in order to design a linear phase high pass filter you need to choose an odd filter length, i.e., a type I linear phase FIR filter.


Defining the SNR or PSNR for color images (3 channel RGB files) In image compression applications, I know that for an $8$-bit $0$ to $255$ level grayscale image the PSNR (peak signal to noise ratio) can be defined as:
$$
\text{PSNR} = 10 \log_{10} \frac{\sum_{m,n} 255^2}{\sum_{m, n}|s(m,n)-\hat{s}(m,n)|^2}
$$
where $s$ is the original 2D image and $\hat{s}$ is the reconstructed image (decompressed). My question is how does one define PSNR (or even SNR for that matter) in 3 color space like a 3 channel RGB image, in which each pixel is associated with a 3-tuple $(R, G, B)$, where $R$ is a Byte representing the value of Red channel and so forth, assuming that the image is an 8-bit depth image.

Your PSNR formula has in the denominator a sum of squared distances in a color intensity 1-D space between original and reconstructed pixel values, the sum is taken over all the image pixels. Therefore, you question can be re-formulated as "how one defines a distance in a color 3-D space".
You did give an implicit answer, one of many possible and still valid, in your question: your elaboration on the definition of color space ("each pixel is associated with a 3-tuple (R,G,B)") suggests that you consider the RGB color space. The Euclidian metric can be selected and you have the formula for a color distance between two color vectors F1 and F2
$$
Δ_{RGB}(F_1, F_2) = \sqrt{(R_1-R_2)^2+(G_1-G_2)^2+(B_1-B_2)^2}
$$
and it is true, but not the whole truth.
The practical usefulness of the color distance definition is task specific. From one point of view, it depends on what is the purpose you assign to your image analysis instrumentation. Considering the image impairment by a reconstruction noise, you may want to select a color space adapted to human vision, a perception based color space, as a hue-saturation-intensity (HSI) color space. The metric in HSI used for color distance calculations is
$$
Δ_{RGB}(F_1, F_2) = \sqrt{Δ_I^2+Δ_C^2} \\
Δ_I = |I_1-I_2|; Δ_C=\sqrt{S_1^2+S_2^2-2S_1S_2cosθ} \\
θ = (|H_1-H_2|<π)?(|H_1-H_2|):(2π-|H_1-H_2|)
$$
In the other project, you may want to construct an image from measurement data for presentation in pseudo colors. It will require to work with the proprietary color space that reflects transfer functions of your detectors.
Summing up, your project may require a standard color space, and you simply select the best color space for your task, or it is a sophisticated research work, and to define a PSNR formula for your image processing is a work of its own.
About "color spaces": a search for "distance in color image space" gives as a first entry an article that I readily recommend to you: https://profs.info.uaic.ro/~ancai/DIP/laborator/lab7/color%20spaces%20&%20distances.pdf
In the article, a section 3.4.3 Opponent Color Spaces (page 62, or 26/34 in the pdf) may be what is of special interest for you.


is root mean square possible to calculate similarity of FFT result? I need to compare FFT result from audio files,
I have 2 FFT result from 2 audio files..

FFT 1 (sample)

-0.16314493488504767 - 0.103707391105263i
0.07863935536550609 - 0.30111206509352917i
0.06753127619272284 + 0.02837438103569126i

FFT 2 (sample)

-0.006953384463464084 + 0.0291883094944081i
-0.09447919915711694 - 0.22676541801892i
-0.23493813662827812 - 0.07765408995141115i

Then I need to calculate the deviation of each real-part point (without imaginer part) from the FFT result then make the result absolute, something like this

-0.16314493488504767 - (-0.006953384463464084)
--> then use Math.abs to get absolute value

then sum all the deviation point, and use Root Mean Square formula to follow with to get the difference value...

where N is my numbeofFFTPoint and the following x is the sum value of my deviation result.

But then, someone told me to use Zero crossing rate formula to get more accurate result

Honestly I don't really understand about Zero Crossing Rate,

what have I tried :

I compare the exact same 2 audio files, then get the deviation and sum result --> 0.0  and of course the Root Mean Square Formula result is 0.0 too..
since there's no difference between the files the it would return 0.0 right ??

question

is it possible just to user Root Mean Square ??
Any advice which formula should I use for this case ?

Thanks :) (sorry 4 bad explaining)

The Discrete Fourier Transform is a unitary transform (maybe scaled).  As such, it preserves the Frobenius norm of the signal.  That is essentially the same as the Root Mean Square calculation.

In other words:

 rms( x - y ) == C * norm( fft(x) - fft(y) )

for an appropriate constant C depending on your fft implementation.

In Matlab or octave, C=1/nfft:

 n=8192;
 x=randn(n,2)*[1;j];
 y=randn(n,2)*[1;j];
 sqrt( mean( abs(x-y).^2 ) ) - norm(fft(x) - fft(y)) / n

 ans = -5.3291e-15


What is the uHz-rotator algorithm? According to this page, the Baudline software tool implements

  ...the uHz-rotator algorithm that operates in the complex domain to
  calculate an extremely accurate estimate of the frequency.

Furthermore

  Under good conditions (strong signal and low noise) the frequency
  accuracy can exceed 8 significant digits (~.000001 Hz).

I tested it with a signal generator (ToneGen). For 82.41 Hz (E2) it displays 82.40999 with averaging and 82.4 without (the rest of the digits keep spinning). The FFT size is only 2048 and the sample rate is the standard 44100Hz.

To compare, an implementation I'm working on, which employs quadratic interpolation, has a poor estimate of anywhere between 84 and 86Hz. The phase difference estimator isn't doing much better. In both cases I only get the precision to the first decimal when I increase the FFT size to 8192.

What is this uHz algorithm and how do you get such precision?

From the (limited) description the uHz rotator algorithm sounds like one of the phase-weighted averages from this site, but it's not an algorithm I am familiar with.

The Cramér–Rao lower bound$^1$ for estimating the frequency of sinusoid with amplitude $A$ in white noise with variance $\sigma^2$ is given by:

$$
\mathrm{var}(\hat{f}) \ge \frac{12}{(2\pi)^2\eta N(N^2-1)}
$$

where $N$ is the number of samples and $\eta = A^2/(2\sigma^2)$ (taken from Kay's Fundamentals of Statistical Signal Processing, Estimation Theory equation 3.41).

So, assuming $A =1$, $\sigma^2 = 0.1$, and assuming that frequency accuracy can exceed 8 significant digits means that $\mathrm{var}(\hat{f}) < (10^{-8})^2$ we get:
$$
\begin{array}
\ (10^{-8})^2 &\ge& \frac{12}{(2\pi)^2 5 N(N^2-1)}\\
\frac{20}{12} \pi^2 \times 10^{-16} &\ge& \frac{1}{N(N^2-1)}\\
N(N^2-1) &\ge& 6.079271 \times 10^{14}\\
N &\ge&  84713
\end{array}
$$

This means that your data length would need to be 84713 samples to achieve this accuracy.

For different values of $\sigma^2$ see the following table.

$$
\begin{array}
\ \sigma^2 & \eta & 10\log_{10}(\eta) & N\\
\\
0.1&    5.\ \ &    6.9897\ \ &    84714.\\
\\
    0.01&    50.\ \ &    16.9897\ \ &    39321.\\
\\
    0.001&    500.\ \ &    26.9897\ \ &    18252.\\
\\
    0.0001&    5000.\ \ &    36.9897\ \ &    8472.\\
\end{array}
$$

PS: I have emailed the author of that web-site and asked for their comment. Let's see if they respond.

$^1$ The Cramér–Rao lower bound is the lowest achievable variance of an unbiased estimator.


Does the impulse response of a system have any physical meaning? In other words, what does the impulse response tell us about the characteristics of the system? For simplicity, let's assume we deconvolve a discrete output and a discrete input to obtain an impulse response that has n points. Does the amplitude of the first (or second, or third...) point indicate something? Does the length n tell us anything?

The impulse response of an LTI (linear time invariant) system has the same physical meaning of the differential equation of a general dynamical system has.

So every output can be computed from the input and the impulse response convolution, without having to solve the differential equation. The frequency charactheristics can also be deduced from the Fourier transform of the impulse response.

Several properties of the (LTI) systems can be seen in the impulse reponse (and Fourier transform of it) expressions, such as slow system, fast system, wide band narrow band systems, overshoot, oscillations etc...


How to calculate the Diagonal loading factor evaluate calculate the inversion of a covariance matrix I am programming a Generalised Likelihood Ratio Test (GLRT) detector. When it comes to inverting a covariance matrix $Ri$, I need  to do a diagonal loading to fix the problem of sigularity of this matrix, but I don't have any info about the technique or algorithm that allows to estimate the right diagonal loading factor factor $\sigma$.

Rather than calculate inverse ($Ri$) I have to calculate inverse ($Ri  +\sigma \cdot I$) where $Ri$ is the covariance matrix.

Need help, please

Short answer: just use $\sigma = 10^{-8}$.
Covariance matrices have eigenvalues $\geq 0$ (theoretically),
so $Ri + 10^{-8} \, I$ will have eigenvalues $\geq 10^{-8}$,
safely non-singular.

A longer answer:
split Covar = S + N, "signal" + "noise", by eigenvalues or by SVD,
Singular-value decomposition
aka PCA, Principal component analysis.
This has several advantages:

gives you an idea of how noisy your Covariance is -- a good thing to know
you can choose a cutoff, e.g. 98 % of the sum of the eigenvalues,
or a number of terms, e.g. 3, to plot and look at
fast, rock-solid implementations:
in Python,
numpy eigvalsh,
numpy svd
and sparse
scipy svds
(sparse, dense or LinOp).

Be careful with large covariance matrices:
eigenvalues can be < 0 because of roundoff,
and they can be very noisy -- see
here
on stats.stack .

See Press et al. for a good discussion of SVD:
"Solution by use of Singular-value decomposition",
Numerical Recipes p. 793 ff.


Can you turn a square wave into a sine wave using a low pass filter? And if it could, would it make the sound of the square wave thinner than before because of losing its harmonics?

To be able to analyze what a low pass filter does first you would need to understand what a Fourier transform is, hence some theory first.
The Fourier transform essentially represents the time-domain information stored in some function (a square wave in your case) in terms of frequencies.
A simple example would be a sine wave $sin(2\pi f_{c} t)$ which contains just a single frequency, $f_{c}$. In the frequency domain, after taking the Fourier transform, it produces 2 frequencies at $f_{c}$ and $-f_{c}$ (the negative value occurs because of the mathematical representation of the Fourier transform using complex exponentials, but nevertheless the positive frequency $f_c$ is what one would expect).
Sine wave frequency response http://www.gaussianwaves.com/gaussianwaves/wp-content/uploads/2013/12/Sine_wave_power_calculation_matlab_2.png?02772a
If you are clear with that, then let's move on to what a filter is. A filter modifies the input and generates some output from it and is often represented in the frequency domain => A filter modifies the Fourier transform (aka frequency spectrum) of an input. Basically, it may attenuate some parts of the frequency spectrum and magnify other parts, and may add a "phase change" (frequencies manifest themselves as complex numbers, and phase changes become relevant then)
Now, lets see what exactly a low pass filter does to a square wave.
A low pass filter can be understood as a filter that cuts-off frequencies above a certain threshold when applied to some input. High frequencies which correspond to rapid transformations in the time domain are removed from the input and only the lower frequencies remain.

(source: www.sfu.ca)
In the image attached above, the function on the left is a square wave and the one on the right is the magnitude plot of its corresponding Fourier transform (as mentioned before, the frequencies are actually complex numbers and the entire Fourier transform would be complete only with the addition of a "phase plot" as well).
Observing the magnitude plot, if one were to pick a low pass filter with cut-off frequency between that corresponding to bin #1 (1st harmonic: 100 Hz) and bin #2 (3th harmonic: 500 Hz), the other higher frequencies would get cut out and there would only be 2 frequencies that remain: -100 Hz and 100 Hz which would correspondingly represent a sine wave (as mentioned above).
Answer:
It is therefore viable to use a low pass filter to create a sine wave out of a square wave. However one is restricted to a sine wave with the frequency the same as the period ($1/f_o$) of the input square wave. Higher harmonics cannot be generated (a band-pass filter would be required for the same).
Image sources: google
EDIT: Thanks to @Matt for pointing out that the question says square wave and not rectangular function. Edited answer to reflect the same.


Nyquist Maximum data rate formula for PCM The Nyquist maximum data rate formula for a binary PCM is given by

$$C = 2B\log_2L$$

I'm not very sure what "$B$" is here. Is it the bandwidth of the signal being sampled or is it the channel bandwidth? Also, is "$L$" the no. of quantization levels or is it the no. of M-ary symbols used for baseband transmission?

I'm also not sure about what it means when we refer to "Bandwidth of a PCM system"? Are we talking about the channel bandwidth or the bandwidth of the PCM signal itself? I looked up these things on internet and there are different explanations at different places.

Can someone please clarify this or at least point me to a reliable source where it is explained without any ambiguity! Thanks!

The first document (from mason.gmu.edu) gives a confusing explanation of the formulas. The second document is what's the standard way of referring to the channel capacity

Now the two formulas are:

$C = 2B ~\log_2(M) ~~~~~~~~~~~~,~~~\text{Nyquist}$
$C = B ~\log_2(1 + \text{SNR}) ~~~,~~~\text{Shannon-Hartley}$

Eventhough the first formula, (referred to as Nyquist in the first document), is assumed to yield channel capacity (of a noiseless! channel which is infinite) it's actually giving the necessary minimum data bit-rate to represent an analog message signal of bandwidth $B$ Hz, and quantized to $M$ levels ($\log_2(M)$ bits) using a PCM representation. According to the Nyquist sampling theorem, the minimum allowed sampling rate would be $2B$ Hz, hence there will be $2B$ samples per second each quantized to $M$ levels, yielding a total of $ D = 2B ~ \log_2(M)$ bits per second data rate to represent the analog message signal. Hence it's not channel capacity, it's the minimum necessary bitrate to represent an anlog signal using $M$-bit PCM technique.

The second formula is what we call: maximum channel (reliable) transmission capacity in bits per second rates, or Channel-Capacity in short, over a more realistic noisy channel (such as telephone wires, all sorts of RF channels, and also the fiber optic or satellite links by modeling them with the inherent noise floors). It gives you the maximum achievable rate of information (or data) transmission through a channel of bandwith $B$ Hertz, by using a modulation scheme with the given SNR (ratio of transmitted signal power to the channel noise power).

Note that despite being treated in communications sections, PCM is not a data transmission format but a data representation format. Analog data is converted to digital by samlping and quantizaiton through an ADC of $N$-bit quantization at a sampling rate of $F_s$ samples per second. The resulting sampled-quantized data is practically represented in binary digital system for processing and transmission. In such a case every $N$-bit quantized sample of the message signal is laid out as $N$ consecutive binary digits (bits) to represent a single PCM sample of message data. When it comes to transmission, those bits can be transmitted one by one as is done by using BPSK or in groups of M-bits as is done by 8-QAM (where 8 PCM bits are transmitted at once by a single 8-QAM transmission symbol).


If we take Z-transform of a signal & find its ROC. How to use this ROC? There are ∞ signals in ROC, suppose we choose any signal what to do with it? If we take Z-transform of a signal & find its ROC. How to use this ROC? There are infinite signals in ROC, suppose we choose any signal from ROC, what to do with it?

Of particular interest is the unit circle ($z=e^{j\omega}$) when it falls within the region of convergence, since the Z-transform $H(z)$ with $z=e^{j\omega}$ is the Fourier Transform. If we compute the Z-Transform of the system’s impulse response (unit sample response), this will give us the frequency response of the system. We can thus use the ROC to determine if the Fourier Transform even exists.
The Z-transform of a discrete function is only unique if the ROC is also provided. Thus, and importantly, by knowing the ROC we can uniquely determine the original discrete function by using the inverse Z-transform. If the system is known to be causal then the ROC is known to be the outside of a circle which resides on the inner-most pole in the z-plane. This property of the relationship of the ROC allows us to mathematically work with non-causal signals by specifying the ROC as a ring.
The actual locations of the singularities (poles and zeros) provides a lot of additional information about the system not provided with the Fourier Transform alone, and these singularities would be the primary locations of interest within the ROC.
Here is another perhaps more intuitive view that expands on the unit circle of the Z-transform being the Fourier Transform and the concept of convergence with regards to the Z-transform result: Every circle centered on the origin within the ROC is the result of the Fourier Transform of the original function multiplied by a real decaying or growing function given by $|z|^{-n}$ where $n$ is a sample number. Thus if the unit circle is not in the ROC, we can determine from ROC a weighting function that we can use in order to determine the Fourier Transform. We see this by comparing the formulas for the z-transform with the Discrete-Time Fourier Transform (DTFT) for an arbitrary function $x[n]$:
DTFT
$$X(\omega) = \sum_{n=-\infty}^{\infty}x[n]e^{-j\omega n}$$
Z-Transform
$$X(z) = \sum_{n=0}^{\infty}x[n]z^{-n}$$
$z$ is a complex variable that has a magnitude and phase given as $|z|e^{j\omega}$ (in general the form $Ae^{j\phi}$ is $A\angle \phi$ for real $A$ and real angle $\phi$), and thus the z-transform is related to the DTFT as:
$$X(z) = \sum_{n=0}^{\infty}x[n]z^{-n} = \sum_{n=-\infty}^{\infty}u(t)x[n]|z|^{-n}e^{-j\omega n}$$
Where $u(t)$ is the unit step function.
Where we see that the z-transform is the Discrete Time Fourier Transform of the causal function $u(t)x[n]$ weighted by the real function $|z|^{-n}$ which is  decaying or growing with $n$ (or constant if on the unit circle).
Following an example provided by the OP in the comments:

Consider a signal having Z-transform X(z) = 3+(2/z) & its ROC is whole
z-plane except 0. Suppose z=1, so X(z) =5

With the explanations above we can see intuitively what is occurring and what this means. First, consider how the z-transform on the unit circle is a correlation of our original function $x[n]$ with the set of complex rotating phasors (individual frequency tones) given by $z^n = e^{j\omega n}$, given that correlation is simply a sum of complex conjugate products. When $z=1$ we are correlating our function $x(t)$ with $z^{n}= 1$ for all $n$, basically a DC signal. So $5$ is the DC response of our system, assuming $x[n]$ was the impulse response (unit sample response) of the system, as then the FT is the frequency response and this would imply that a DC signal at the input to the system would have a gain of $5$. If $x[n]$ wasn't an impulse response but an arbitrary signal, $X(z)=5$ for $z=1$ would imply that the DC content of the signal $x[n]$ (which is the mean value) would be $5$ in proportion to the levels everywhere else for $z$ on the unit circle (basically the FT at $f=0$ or DC, which is a correlation of our signal to all possible frequencies).


MIT exercise 6.003 HW2 - Concept of system initially at rest I am following the MIT open course you can find here. My question is about one of the exercises given as homework in the latter and more specifically I think I am missing something on the concept of "system initially at rest".

Homework [exercise 8]: HW2
Solutions: SOL2

In the proposed solution he is solving the problem using differential equations. What I am trying to do is to use functionals (as this is new to me). I think I end up missing an integration constant due to the hypothesis of "system initially at rest". My solution goes as follows:

System differential equation:
$$ \ddot{y}(t) = x(t)-y(t)$$

We seek the response to a step function i.e. $x(t) = u(t)$ where $u(t)$ is zero for $t<0$ and one elsewhere. Now the system functional is:

$$\frac{Y}{X} = \frac{A^2}{1+A^2} = \frac{i}{2}\left( \frac{A}{1+iA}+ \frac{A}{1-iA} \right)$$

where $Af(t) = \int_{-\infty}^{t} f(\tau) d\tau $ is the accumulator operator as defined in the course. Given the impulse function(al) $\delta(t)$, knowing that $u(t) = A\delta(t)$ and using the commutative property of operators:

$$\frac{Y}{X}u(t) = \frac{iA}{2}\left( \frac{A}{1+iA}+ \frac{A}{1-iA} \right)\delta(t)$$

The fundamental mode being $e^{pt}$ where $p$ is a pole:

$$\frac{Y}{X}u(t) = \frac{iA}{2} \left( e^{-it} - e^{it} \right)u(t) = \frac{i}{2}\left( \frac{e^{-it}}{-i} -  \frac{e^{it}}{i}\right) = -cos(t)u(t)$$

While the correct solution is (physically obvious): $\left(1-cos(t)\right) u(t)$. I don't understand where I am doing wrong except for the initial hypothesis "system at rest". Can you tell me where the mistake is and, if the hypothesis is wrong, tell me exactly in what conditions it holds?

The mistake is in the final integration of the term

$$\frac{1}{2i}\big(e^{it}-e^{-it}\big)u(t)=\sin(t)u(t)\tag{1}$$

because you seem to ignore the step function $u(t)$. Integrating $(1)$ gives

$$\int_{-\infty}^t\sin(\tau)u(\tau)d\tau=u(t)\int_0^t\sin(\tau)d\tau=\big[1-\cos(t)\big]u(t)\tag{2}$$


what is the correct order of applying a nonlinear operation and downsampling (with anti-aliasing) In UltraSound imaging, after beam-forming and obtaining the envelope, Log-compression which is a non-linear operation and scan conversion which is a 2D interpolation (with inherent down-sampling) is performed. I was wondering what would be the correct order of log-compression as a nonlinear operation and the 2D interpolation which inherently down-sample the signal.
Generally if we perform some non-linear operation and down sampling what would be the correct order?

When one nonlinear and one linear operation are suitably adapted to the problem, the nonlinear one is quite often applied first. Nonlinear processing is often applied to modify data such that more classical algorithms apply more easily, even when the original data does not follow the proper assumptions. Examples are:

removing outliers/trimming data before computing averages,
applying variance stabilizing transformations (Anscombe) to use regression needing constant variance,
computing a logarithm to convert a multiplicative problem before using an additive technique on the log-transformed data.

Of course, applying a non-adapted nonlinear transform first can be very harmful. Yet in you case, a log-coompression may either:

reduce the impact of extreme data values
stabilize an heteroskedastic noise
account for multiplicative effects

before an interpolation (a kind of weighting average or filtering). This goes under the hood of homomorphic filtering.


How to perform carrier phase recovery in software? What are the options for recovering the phase of a BPSK signal in software? The only resources I can find online contain circuit diagrams – there seems to be no one willing to explain this to people unfamiliar with analog circuits.

I'd like any explanation of how carrier recovery works in theory, and would also love pseudocode or code examples.

To demodulate a phase-shift keyed signal, of which BPSK is the simplest, you have to recover the carrier frequency, phase, and symbol timing.

Bursty Signals
Some signals are bursty and provide a known data sequence called a preamble or mid-amble (depending on whether it shows up at the beginning or middle of the burst).  Demodulators can use a matched filter that "looks" for the known data sequence and use it to figure out the burst's frequency, phase, and symbol timing.  They do this for every burst and don't need to bother with "drift" (the gradual distance that builds up between the demodulator and the received signal as small errors accumulate) because the bursts are generally short enough that the drift isn't a problem.

Continuous Signals
Then there are continuous signals.  They are in many ways tougher to recover than bursty signals because you normally do not have a known data sequence to aid in locking up with the signal, and you do have to worry about drift even after you lock onto the signal.  I'll try to describe at a high level the main steps that are typically used to recover continuous signals.

Carrier Recovery

Usually you know what frequency the signal you are looking for will be at, or at least what set of frequencies it could be at.  Even with this knowledge, though, you usually have to be able to correct frequency offset because no two transmitters transmit at the same frequency.  There is always some error.  The usual method, then, is to mix the frequency that you think that the signal will be at, and then correct the residual frequency error.  This can be done with a Costas Loop, or by taking the fourth power of the baseband signal data and looking for a frequency spike.  There should be a frequency spike at the carrier offset * 4 (e.g. if you take an FFT of the data to the 4th power and see a frequency spike at 8300 Hz, that means the carrier offset is 8300/4 = 2075 Hz).  This is a very effective means of getting an initial lock on the frequency offset.  You can also use it to compensate for drift if you re-do it occasionally.  There is another way to compensate for drift which I will touch upon later.

Carrier Phase

At this point if you plotted your complex data in the complex plane (x-axis is real, y-axis is imaginary) it should look something like the following-

If you look closely you can see two dense areas towards the ends of the blurry line.  Those are the BPSK constellation points.  The points in between are the transitions between the constellation points.  Those will clear up once we get the symbol timing.  The reason that the line is at an angle is because of the carrier phase.  This can be measured by mirroring the signal by multiplying all points that have negative real values by $e^{j*\pi}$

and then taking the mean of the angles of the points.  Once you have calculated that, subtract that angle from all points that have had the carrier offset removed by multiplying the points by $e^{j*-\omega}$.  You can also compensate for carrier offset drift with this technique by dynamically updating the phase offset.  Once the data is phase corrected it should look something like this-

Once the data is phase corrected you can drop the imaginary part of the data because it does not add any information.

Symbol Timing

You should normally know, a priori, the symbol period of the signal you are trying to demodulate.  If you do need to determine the symbol period/frequency though, you can do it in a similar way to how the carrier offset was detected.  You can square the data which will cause a frequency spike at twice the symbol frequency.

Like with the carrier you will have to get the phase (timing) right and then compensate for drift.  The usual method for both of these problems is to look for the zero crossings.  Unless the noise is quite bad it should only cross the zero-point in the middle of a symbol transition from -1 to 1 or from 1 to -1.  Even if noise does cause this to happen in the middle of a symbol it won't happen very often.

The picture above is usually called an "eye diagram" or "eye pattern".  It is two symbol periods wide, and has many symbols "stacked" on top of one another.  I don't know if you are familiar with oscilloscopes or not, but you can get an oscilloscope to show a picture like this.  Anyway, the two "X's" are symbol transitions.  The high to low lines in the X's are when the symbol transitions from a 1 to a -1, and the low to high lines in the X's are when the symbol transitions from a -1 to a 1.  The point in between, where the pink line is, is the optimum place to sample the data to see if the symbol is a 1 or -1.

This same technique can be used to handle symbol timing drift.  Take a running average of the distance from the prior zero crossing and following zero crossing.  If the two averages are about the same, all is well.  If one is bigger than the other then you need to shift where you are grabbing your sample.

Once you have sampled the symbols at the correct points your data points should get something that looks like the classic BPSK constellation points.

I hope this helps.


How should I select window size and overlap? I'm new to DSP and I have a wav file with a number of different frequencies playing at different volumes. I'd like to plot the power of the 350Hz frequency signal at each point in time. The source is rotating so should see the power rising and falling periodically. To start with I broke the signal up into windows and then implemented Geortzel's algorithm (from Wikipedia). This seems to work: https://www.dropbox.com/s/aq1aigdx8dq306q/Screenshot%202015-03-18%2017.38.20.png?dl=0 but for different window length and window overlaps I get different plots.

I'd love some guidance on what the correct way to choose the window length and the amount by which each window overlap is? window_length and window_skip in the code.

Also should the windows be processed (i.e. convolved with a window function like Hanning) I've read about this but I'm not clear on what this would do to the output, nor how to choose the length of the window function?

Here's the audio file: https://www.dropbox.com/s/knvwqwg7s3mcz56/example.wav?dl=0

Here's my code, comments welcome:

import pylab
import scipy.io.wavfile
import numpy as np
import math

def goetrzel(data, target_frequency):
        s_prev = 0
        s_prev2 = 0
        normalized_frequency = 2.0 * np.pi * target_frequency / len(data)
        coeff = 2.0 * np.cos(normalized_frequency)
        for sample in data:
                s = sample + coeff * s_prev - s_prev2
                s_prev2 = s_prev
                s_prev = s
                power = s_prev2 * s_prev2 + s_prev * s_prev - coeff * s_prev * s_prev2
        return power

def sliding_window(data, length, skip=None):
        if skip is None:
                skip = length
        n = len(data)
        for i in range(0, n - int(length), int(skip)):
                yield data[i:i+length]

detect_frequency = 350.0
#detect_frequency = 485.0
#detect_frequency = 633.0
#detect_frequency = 796.0

rate, data = scipy.io.wavfile.read('example.wav')

window_length = 2 * detect_frequency # How should I pick this?
window_skip = window_length # How should I pick this?

power = [goetrzel(window, detect_frequency) for window in sliding_window(data, window_length, window_skip)]
pylab.figure()
pylab.plot(power)
pylab.show()

Using overlapping rectangular windows has the potential to give an extra 'weighting' to certain sections of the audio unless you choose the overlap carefully. For example, with an overlap of 25% you're accounting for half of the data twice and the other half just once. This is not necessarily an issue depending on the application, but is something to bear in mind. Choosing overlaps of 1/2, 2/3, 3/4... solves this problem (assuming a signal of reasonable length).

For a given overlap, the choice of window length affects two things really - temporal resolution and coefficient quality.

The first is slightly easier to explain and visualise. If the window increases in length then you have fewer windows spread across your audio and therefore a greater time interval between each - i.e. lower temporal resolution.

The coefficient quality, however, increases with window length. If you perform a Fourier transform over a longer sample you get a higher frequency resolution and each coefficient is therefore more 'representative' of its frequency. Think of a coefficient as a frequency bin rather than a single frequency value. Shorter windows leads to wider frequency bins and smoothing/smearing of the coefficients, whereas longer windows leads to narrower, more 'distinct' bins. A quick way of working out the bin size is to do 1 / (length of window in seconds). If this comes out to be 20 Hz and you want to be able to compare the coefficients at 330 Hz with those at 340 Hz then you might consider upping the window length.

You can probably see the tradeoff in window length here - longer windows give less temporal resolution and better coefficient quality; shorter windows give higher temporal resolution and lower coefficient quality.

Currently you're using not using a window function, i.e. your windows are rectangular. Window functions are useful in that they can make your window of data appear more periodic than it actually is. Remember that the Fourier transform assumes periodicity. Any significant differences in the start and end of your window of data become a sharp step change under that assumption, producing a lot of extra coefficients/noise. Multiplying the data by a window function can help mitigate that problem to some extent by reducing the magnitude down to zero or near-zero at the edges of your window. Common window functions used for audio are the Hann and Hamming windows.

What this all means for you:

I'd recommend using a window function to reduce noise if your windows of data are non-periodic, which is likely. Start out with a Hann window.
Id also recommend an overlap of 50%. When used with a Hann window this value has the advantage that 50% overlapping Hann windows sum together to a constant magnitude of unity. An overlap of more than this results in better coefficient quality at the expense of extra processing, but 50% is a good starting point.
With the overlap chosen your decision about window length should be based on the tradeoff mentioned earlier and depends on your application. If you want a high temporal resolution then you're going to have to sacrifice quality. If you want to be able to compare tightly spaced frequencies over time then you'll have to sacrifice temporal resolution.


limit of the integral in the anti-transform if I have an autocorrelation on a limited domain

r($\tau$)=$\lim_{T \to \infty}\frac{1}{T}  \int_{0}^{T} u(t)u(t+\tau)dt$

The power spectral density, obviously will have an infinitive domain:

S(f)=$\frac{1}{2\pi} \int_{-\infty}^{\infty} exp(-i\tau 2 \pi f)r(\tau) d \tau $

If I do the anti-transform of the PSD, which will be the limit of the integral?

from -$\infty$ to $\infty$ or from $0$ to $T$??

Thank

I think you are making this more complicated than it needs to be.

The Fourier and autocorrelation integrals are defined from $-\infty$ to $+\infty$. That's always correct and the saftest way to write them.

If you are integrating over a function with limited support on $[0,T]$ than you can adjust the integration interval to  $[0,T]$ too. The areas outside the support interval is all zero so it makes no difference whether you include them in the integration or not. In this case both versions are identical and both versions are correct. You can use whichever one is more practical for your specific problem or application.

In your specific example: if $u(t)$ has support on $[0,T]$  then it's autocorrelation $r_{uu}(t)$ also has support only on $[0,T]$. That means you can calculate the PSD  by executing the Fourier integral over either $[-\infty,+\infty]$ or $[0,T]$. Both are correct and both will give you the same result.

$S(f)$ has infinite support so you need to use $[-\infty,+\infty]$ for the inverse Fourier integral.

$$S_{uu}(\omega)=\int_{-\infty}^{+\infty}r_{uu}(t) \cdot e^{-j \omega t} d \omega = \int_{0}^{T}r_{uu}(t) \cdot e^{-j \omega t} d \omega$$

$$r_{uu}(t)= \frac{1}{2 \pi} \int_{-\infty}^{+\infty} S_{uu}(\omega) \cdot e^{+j \omega t} d \omega $$


polyphase sample rate conversion with non-integer factor I want to do sample rate conversion by subsequently upsampling with factor I=5, and then downsampling with factor D=9.

I have designed a nyquist sample rate conversion filter h() of length M, with matlab's filterbuilder tool for interpolation factor I=5, and decimation factor D=9. Since D > I, I have chosen the normalized cutoff frequency of the filter to be at 1/9, by setting "Band" under filter specifications to 9.

Then I compute I polyphase filters pk(n), which have length K=M/I by sampling h() according to:

pk(n) = h(k + n*I),   for    k=0,...,I-1    and     n=0,...,K-1.

Then I only compute the output samples for each polyphase filter, which are sampled by the downsampling operation. In other words, I compute the output y[m], by sampling the outputs of the polyphase filters yk[] according to:

y[m] =   yk[ floor( m * D / I  ) ]
k    =   ( m * D ) modulo I

Thus, I do not compute samples of yk[] that are not used in the output.

I apply these polyphase filter's first to do sample rate conversion for the rows of the image. Then I apply the same filters to do sample rate conversion for the columns.

However, I get a distorted output image which still, clearly contains alias:

Does anyone know If I'm conceptually doing anything wrong?

Your approach seems conceptually sound.  I suspect that the problem is one of two things: 1) you have implemented it incorrectly, 2) Your filter bandwidth is too wide.

You are implementing the "smart" approach in that you are doing it the computationally efficient way.  I would try doing it the dumb way (insert I-1 zeroes between every sample, filtering with the entire filter, get rid of the extra samples) to verify your implementation and to see if the filter itself is the problem.  You can verify your implementation by comparing your resulting samples against the samples produced by the dumb approach.  After you have filtered the upsampled data one of the sample "phases" (there are "D" phases whose start samples are samples 0...D-1) should be the same as the smart approach's results.  If not, something is wrong in your implementation of either the smart or dumb approach.

If the results do match, then the problem is your filter.  You need to look closely at the bandwidth to make sure that it isn't too wide.  Since D > I the ideal filter bandwidth, assuming a perfect brick wall filter, is $\frac{I}{2D}$, where the factor of two comes from the Nyquist rate.


Perform inverse distance weighting interpolation using multiple images matlab I have a total of 4 highly identical grayscale images, where image1 is taken as the reference image, the rest of the images are to be used to enhance image1.

My objective is to enlarge the image to 2 times the original size. Instead of taking average values from neighbouring pixels to fill the upsampled pixels, is it possible to do inverse distance weighting interpolation, so that my image remains sharp?

My other concern is that I need to fill the upsampled pixels with the pixel values from the 3 other images. Is it possible to do interpolation with 4 input images?

imresize will enlarge your image, and it gives you a choice of nearest neighbor, bi-linear, or bi-cubic interpolation. As to using the other images for interpolation, you would have to write that yourself.


How does a telephone change your voice? Is it due to vocal length normalization during MFCC? Please be kind enough to explain.

In the case of switched analog POTS phone lines, the audio response is shaped by the limited frequency response of the carbon microphone and the speaker in the handset.


up down sampling and low-pass filters I am trying to figure out what seems very trivial thing about down sampling and up sampling. The problem that I am facing is how to find the effects of up sample on a signal followed by a filter and/or a filter followed by down sampling. Here is a problem related to this from the Oppenhiem and Schafer book.

I am trying to see what happens to x[n] after it pass the hlp[n], considering that hlp[n] has wc = pi/4. Will there be a gain on this or just the frequency is affected?

Please help, Thanks in advance

I think, this will act like low pass filter with $\pi/2$ instead of $\pi/4$ cutoff frequency


What determines peaks in FFT? I ran FFT on three audio files and found that the results for some have more peaks than the other. Could anyone give me any conceptual explanation as to what determines these peaks? Below are plots of the audio sound envelope (left) and the corresponding FFT results (right). If I look at the plots on the right, e.g., the FFT for audio 1 has no peaks between 0.33 and 0.67 but the FFT for audio 2 and 3 do. I am not quite sure whether these additional peaks are meaningful or just artifacts of my particular sound file and FFT setup. More explanations about the audio files below.

Each peak in the audio plots represents a sinewave tone. The smallest inter-tone interval is 187.5 ms. Other inter-tone intervals could be 2x, 4x, or 6x187.5 ms. Each audio file lasts 5x3 seconds. The first row is the result from the exact same 3-second rhythm looped 5 times. The last two are the repetition of different 3-s rhythms. The inter-tone intervals in the three files are as follows:

Here's my FFT script:

n = length(sound);
nyquist = srate/2;
f = linspace(0, nyquist, floor(n/2)+1);
fourier_coef = fft(snd);
fourier_coef = fourier_coef/n;
amp = abs(fourier_coef);
amp = amp .*2;

Thank you.

Teresa

Peaks are determined by energy in a given frequency range which can be more or less visible depending on the dimension of frequency bins.


Matlab - Bode plot of Lag Filter + Integrator I am rather new to Matlab and I just cant make sense of what I see in the bode plot of the continuous and discrete version of the same function. The bode plot of the continuous function looks as expected. However the bode plot of the discrete version has a phase offset of +90 degrees and the gain stays the same at lower frequencies. It's basically a lag compensator with an integrator. This is not the final result I am going for, but the easiest example I could think of to make the problem as clear as possible. Please note that I am using a rather high sampling rate of ~12 MHz, which is what I need in the final hardware design. The high sampling rate seems to be part of the problem. With lower sampling rates the frequency at which the bode plots differ gets shifted to the left.

Can anyone explain this behaviour and maybe how I can achieve a result that is closer to the continuous bode plot using the high sampling rate?

Code of the example:

FS = 48000*256;
T_sample = 1/FS;

accu = tf([0 1], [1 0]);
lag = tf([10 1], [100 1]);

lag_d = c2d(lag, T_sample, "zoh");
accu_d = c2d(accu, T_sample, "zoh");

bode(lag*accu, lag_d*accu_d);

I tried it in Octave, there's definitely a glitch like you said. I tried it with "Tustin" instead of "zoh", same result.

However I was puzzled by your high sampling frequency. Your lag controller has a pole at s = -0.01 and a zero at s = -0.1, they are both "slow" while your sampling frequency is really high.

How to fix it :
Change your sampling frequency to something much lower than 12 MHz, like 10 kHz for instance and replace "zoh" by "tustin". "Zoh" should only be used when converting a "process" transfer function to the discrete domain in order to take in account the effect of the DAC, we don't use it for compensators.

The root cause however is probably some quantization issue in the c2d function.

FS = 10e3;
T_sample = 1/FS;

accu = tf([0 1], [1 0]);
lag = tf([10 1], [100 1]);

accu_d = c2d(accu, T_sample, "tustin");
lag_d = c2d(lag, T_sample, "tustin");

figure
bode(lag*accu, accu_d*lag_d);

Edit 1 :

The issue seems to be with the lag controller, not the accumulator. You should try to move the pole and the zero to the left, this should fix your problem since you must keep your sampling frequency really high:

try this (same DC gain, but pole and zero further to the left) :

lag = tf([1 1000], [10 1000]);


Kalman Filter | Why Shift the Belief Distribution I am looking at the video lecture on Kalman Filter. I got some understanding that if the robot moves that we also shift the enviornment belief model accordingly (the same mean amount the robot moved). I am quite confused that why do we move the belief model of the enviornment. The enviornment itself is stationary so if the robot moves the door will remain in place so their distribution should not be shifted. In the following snapshot, the doors are our measurements and we form the belief of where the robot is with respect to the doors. When the robot moves the doors stay stationary. Why do we move the posterior model of as shown in the image below?

The idea is that the robot is trying to figure out where it is.  The only information it can use is:

Its door sensor which says either you're next to a door or you're not next to a door.
Its map of where the doors are.

That's how we go from the "maximum confusion" (uniform) plot to the one with three humps. This uses the first you're next to a door measurement.

Then, the robot knows it moves to the right.  Therefore it must update its probability with a similar movement to the right.  However, its motion sensors are a little uncertain, so this has the effect of smoothing the three bumps (well, the whole distribution) a little bit. So the three bumps move to the right and they spread out a little.

The peaks of the three spread out humps are still the robot's best guess as to where it is on its map.

Now we get the second you're next to a door measurement.  We multiply our prior (the slightly spread out three-humped distribution) with a new three humped distribution (from our known map).  The prior second-door location and the map second-door location will align, so that will generate the largest peak.


Discrete wavelet transform; how to interpret approximation and detail coefficients? I am trying to understand Wavelet transform. So far I have understood the basic theory of it. But I am not able to get my head around how to interpret both coefficients.

I am using PyWavelets package of python, I have a time-series data for 1 year equally divided into 15-min time interval. I am trying to get high-frequency components for each day. Fourier transform gives the information about high-frequency. But I would like to get time-information also i.e. at what part during the day there occurs a spike(high-frequency). If someone could please help me figure out this and how to plot time-frequency. I have searched a lot but couldn't find a perfect answer. I am not looking for a free solution without trying; Here is everything I have tried

tdayde=d2012_2['Date'].iloc[0]
endde=d2012_2['Date'].iloc[-1]
skiphr=datetime.timedelta(hours=24)
blockde=tdayde+skiphr

I have timestamps attached to data, so to plot for each day I am using

healtharrayde=d2012_2
healtharrayde=healtharrayde[(healtharrayde['Date'] >= tdayde) & (healtharrayde['Date'] < blockde)]
sig=list(healtharrayde['MCP'])
  cA, cD = pywt.dwt(sig, 'db2')
    plt.plot(cD)
    plt.plot(cA)

So if there are 96 values in array, both cD and cA are 48 values each. I would like to understand and plot cD(detail components) with respect to time, So how should I do? because there are 96 timestamps(MCP is column name for values)

Sorry if it's a silly question but I really need help. Thank you

Wavelet transforms can be more difficult to interpret than FFT at face value due to the various representations, nomenclature and output formats. I had to study more than 15 resources to get a good sense of the variety and which one is used by Pywavelets (which does not provide much theory or explanation in its documentation).

In order to grasp the meaning of cD and cA coefficients, it is helpful to run through a basic example wavelet transform calculation.
Here's a simple step-by-step calculation of what happens in a multi-level DWT (your example is basically the first level). In this representation, they concatenate cA and cD coefficients side by side.

When running DWT using a specific basis function, the signal is fed through (inner product with) a high pass filter (difference filter) and a low pass filter (smoothing / averaging filter), each of which is unique to the wavelet basis function.

For 'db2', the high/ low pass filtering each has two terms, and occurs with a step size (stride) of 2, therefore, after the filtering is completed, you also get a downsampling by 2 of the original signal. Actual length will depend on the filter length and the signal extension mode.
The high pass filtered result gives you the cD coefficients.
The low pass filtered result gives you the cA coefficients.

In your example, each of the 48 cD and 48 cA coefficients corresponds to 2 of the original 96 data points.
Therefore, to plot the cA and cD coefficients in time, just reduce time resolution by 2.

Hope this answers your question, if not, here're some more resources:

This has some simple equations and intuitive comparisons to FFT
This has some nice graphical representations and color coding to explain multi-level DWT
Another guide with detailed explanation and illustrations


Calculate the Standard Deviation of Fundamental Frequency (MFCC) I'm implementing a gunshot detector following the article "Algorithm for Gunshot Detection Using
Mel-Frequency Cepstrum Coefficients
(MFCC)" (paywall).

In the article, the authors uses 22 features based on MFCC (coeficients normalized between 0.1 and 0.9). One of the features is "Standard deviation of fundamental frequency". I already spent all my day trying to figure out the proper way to do this calculation.

The expected value for a gunshot is near 0.04.

You calculate fundamental frequency of the signal with some pitch extraction algorithm, for example YIN, https://github.com/ronggong/pypYIN. It is a vector of numbers, then you calculate standard deviation of it, that would be a number.

Overall, this feature is not very important for detection. Also, the paper you are trying to reproduce is not quite modern, these days you can get much better results with VGGish and Audioset, something like

https://github.com/luuil/Tensorflow-Audio-Classification

or

https://github.com/IBM/MAX-Audio-Embedding-Generator


Discrete Hilbert transformation of a spatial signal first of all, sorry for this basic question in signal processing, I'm pretty new in this field and have been stuck in a situation a little beyond my expertise.

The idea is to derive the instantaneous wavenumber direction as an estimation of the direction normal of a 2D field by using a directional Hilbert transform in one of the directions.

As far as I can recall correctly, all I have to do is to do a 1D Fourier transform of each data vector, apply a $\pi/2$ phase-shift in the wavenumber domain and then apply an inverse Fourier transformation to get the Hilbert transform.

So the Hilbert transform is defined as:
$$
H(x(t)) = \frac{1}{\pi} PV \int_{-\infty}^{\infty}\frac{x(\tau)}{\tau-t}d\tau
$$
where $PV$ denotes the principal value of the integral and $H$ is the Hilbert transform of a vector $x(t)$. This is supposed to be equal to the convolution $\frac{1}{\pi t} * x(t)$ and

$$
F(H(x(t))) = F(\frac{1}{\pi t}) F(x(t)).
$$

So the Fourier transformation and multiplication of the two factors on the right side is equal to the Fourier transformation of the Hilbert transformation of my vector, right?

My question now is how to apply the phase shift on the Hilbert transform? is it simply multiplying the product on the right side by $\pi/2$ in the wavenumber domain?

Thanks a lot for any help!
Cheers!

UPDATE:

So I coded the Hilbert transformation in $C$ but unfortunately I do not get a phase shift of $90^\circ$ but somewhat different from that and with a wrong amplitude. My code looks like this:

    void hilbertfun(float *vin, float *vout, int n, int d)
    {
static complex *A;           //working array
int i;
float pi = 4. * atan(1.);
int signfun;
static float *rkx;

void rk(float *ak, int n, float d, int ind);

A  = (complex*) alloc1(n,sizeof(complex));
rkx = (float*) alloc1(n,sizeof(float));

//copy input array to complex array
for(i=0; i<n; i++){
    A[i].r = vin[i];
    A[i].i = 0.;
}

// perform inverse fft
fouriertransform(1,n,A);

// calculate wavenumbers
wavenumbers(rkx,n,d,1);

// check wether wavenumbers are positive or negative >> apply phase shift by using the impulse response of the Hilbert transformer
// + apply phase shift
for(i=0;i<n;i++){
  if (rkx[i] > 0) signfun = 1;
  else if (rkx[i] == 0) signfun = 0;
  else if (rkx[i] < 0) signfun = -1;
    A[i].r = - A[i].i * signfun;
    A[i].i = A[i].r * signfun;
}

// perform forward fft
fouriertransform(-1,n,A);

for(i=0;i<n;i++) vout[i]=A[i].r/n; //write into output and normalise

}

The Fourier transform and the normalisation work perfectly fine without the phase shift (Original (red), Output (green)):

Yet, with the phase shift I get this result (Original (red), Hilbert  transformed Output (green)):

What am I missing?

You need the Fourier transform of the impulse response of the Hilbert transformer

$$h(t)=\frac{1}{\pi t}\Longleftrightarrow H(j\omega)=-j\textrm{ sign}(\omega)$$

So all positive frequencies are multiplied by $-j$ and all negative frequencies are multiplied by $j$. This corresponds to a phase shift of $-\pi/2$.

EDIT:

I'm referring to your updated question including the code. I can't see if the computation of the wavenumber is OK, but one problem I definitely see is in applying the phase shift. You change A[i].r and then you re-use it (as if it were the orignal value) to compute A[i].i. Furthermore, I think the signs are wrong. If signfun=1, A[i].r=A[i].i and similarly for the computation of A[i].i. But the most important thing is to use a temporary variably to store the old value of A[i].r for computing A[i].i.