Why use parametric based estimation methods - confusion regarding terms Using the probability density function (pdf) we can estimate an unknown parameter using methods such as Maximum Likelihood estimation. If the pdf is not available, then Least Squares can  be used. Other methods for pdf unavailable situations are Least Mean Squares.

Then there are other methods which are called as heuristic based approaches such as particle swarm optimization, genetic algorithm, ant colony optimization. These are also applied to find an unknown parameter without using the pdf. This is quite a long post, so  a big thank you and I really appreciate the patience in taking the time out to read and understand my concerns. My questions are:

1) What is the difference between parameter estimation and optimization?

2) Is there a rule of thumb when to use estimation theory and optimization?

Many research articles and text books use these two terms interchangeably based on my understanding. For example, the abstract in the paper titled, "Parameter estimation with bio-inspired meta-heuristic optimization: modeling the dynamics of endocytosis" which is given below:

We address the task of parameter estimation in models of the dynamics of biological systems based on ordinary differential equations (ODEs) from measured data, where the models are typically non-linear and have many parameters, the measurements are imperfect due to noise, and the studied system can often be only partially observed. A representative task is to estimate the parameters in a model of the dynamics of endocytosis, i.e., endosome maturation, reflected in a cut-out switch transition between the Rab5 and Rab7 domain protein concentrations, from experimental measurements of these concentrations. The general parameter estimation task and the specific instance considered here are challenging optimization problems, calling for the use of advanced meta-heuristic optimization methods, such as evolutionary or swarm-based methods.

However, in signal processing, optimization is rarely combined or taught in conjunction with parameter estimation. For example, in the book by Steven Kay titled, "Fundamentals of Statistical Signal Processing, Volume I: Estimation Theory" never mentions optimization as another method to do estimation. However, Gradient Descent is an optimization algorithm which uses the derivative of the function and equates to zero. If this is the usual way of doing optimization, then Maximum Likelihood estimation also find the derivatives with respect to the unknown parameter using the log-likelihood function, which are equated to zero to obtain the estimates. This brings to my last 2 questions which are:

3) Is MLE an optimization or estimation technique?

4) Heuristic approaches such as particle swarm etc do not find derivatives to obtain the minima or maxima of the objective function. Then why are they called optimization methods if they are not optimizing the cost function in the usual way which is taking the derivatives and equating them to zero?

Hi: I'll try to answer as briefly as possible and only with respect to statistics. not dsp.

In statistics, if you have a nice pdf such as the normal distribution, then maximizing the likelihood is equivalent to minimizing the sum of squares of the residuals ( often called errors ).

In other cases, where you either have a complicated distribution ( maybe the model is non-linear and results in a more complex likelihood), then maximizing the likelihood is not equivalent to minimizing the residual sums of squares. So, you may still choose to maximize the likelihood but because the likelihood is complex, one often needs to turn to a numerical optimization algorithms such as gradient descent, BFGS.

So, I think your confusion is due to the fact that it ( maybe unfortunately ) happens to be the case that maximizing the likelihood of the normal is equivalent to minimizing the residual sums of squares.

Then there are other cases maximizing the likelihood is too difficult or the likelihood is not even derivable in closed form.  In this case, people define other objective functions ( heuristics ) that still measure how "well" the model fits the data. Doing that, allows one to use the genetic, particle swarm type optimization schemes that you refer to. These schemes are often called non-parameteric because these objective functions are not tied to a parametric likelihood ( such as the normal,t, gamma etc ). In fact, the likelihood is often not even known.

Basically, what type of optimization one uses depends on whether you have a pdf and how complicated it is.

A reallly rough categorization is as follows:

1) normal case: : straightforward: maximize it or minimize residuals squared.

2) non-normal case: likelihood is complex, maybe non-linear but known  : bfgs or gradient descent or some other numerical method that can handle the
closed form version of the likelihood.

3) likelihood too complex or unknown or doesn't exist: define one's own objective function that measures the quality of the model and use algorithms such as particle swarm etc.

Note that in 2) above, you have to worry about whether the optimum you hit is local, global etc and this can be extremely tricky.
This may be the case in 3) also. I'm not as familar with genetic algs as the numerical methods.

I hope this helps some. I'm not sure of a good reference for what you asked about. Kenneth Lange has a book "numerical methods in statistics' that may be relevant but I don't have it. Good luck.


How much of a problem is the Gibbs phenomenon and would this solve it? By iterating or nesting a sine curve like in this question

I get curves like these:

that seem to tend to a square wave.

The eight case of these looks like:

which here on purpose was chosen for it's roundness.

The Mathematica code for these plots can be found here:
http://pastebin.com/6UK1u1uX

I don't know much about signal processing but I recalled the Gibbs phenomenon in square waves after I saw these curves.

  Would they solve the problem with the Gibbs phenomenon in the case of square waves?

In the Fourier transform this kind of function is not of any use though I understand.

Edit 13.1.2013:

Sawtooth wave:
http://pastebin.com/JNg7bzzB

Triangular wave (partial sums instead of integrals):
http://pastebin.com/wRCBV7NF

Dirac comb http://pastebin.com/QMSMQf26

I would say that this is interesting.  There has been a lot of work done regarding the study of Gibbs Phenomenon.  You should check out the following document to get a better understanding of how it comes up in practical DSP applications:

http://people.clarkson.edu/~ajerri/books/examples/Gibbs_Book.pdf

The typical way to manage Gibbs Phenomenon is to use time domain window functions that taper data at the start and finish.  The window functions reduce spurious contributions to frequency domain information that come from discontinuities at the edges of data sequences.

I haven't seen much application of generating signals by composing them with individual sine waves.  Generally signal generation is done directly in the time domain.  I'm not sure how the function constructions you've documented can be employed to solve a practical problem, but perhaps there is application if you can identify an appropriate problem.


Discrepancy when calculating LTI system output using inverse z-Transform I'm given a difference equation, $y[n]-0.4y[n-1]=x[n]$, and asked to find the natural response $y_n[n]$, forced response $y_f[n]$ and complete response $y[n]$ if $x[n]=4 (0.25)^nu[n]$ and $y[0]=0$.

By one approach I've read,
$$
\text{characteristic equation:}\quad z-0.4=0 \quad\therefore y_n[n]=A(0.4)^n \\
\\
y_f[n]=B(0.25)^n \quad\text{substitute back into diffeq:}\\
B(0.25)^n -0.4B(0.25)^{n-1}=B(0.25)^n - \frac{0.4}{0.25}B(0.25)^n= 4(0.25)^nu[n]\\
\therefore B=-\frac{20}{3}, \quad y_f[n]=-\frac{20}{3}(0.25)^nu[n]\\
\therefore y[n] = A(0.4)^n-\frac{20}{3}(0.25)^nu[n]\\
y[0]=0 \rightarrow A=\frac{20}{3} \\
y[n]= \frac{20}{3}(0.4)^n-\frac{20}{3}(0.25)^nu[n]\\
$$
Great! Alternatively, I should be able to determine $y[n]$ by taking the inverse z Transform of $Y(z)=H(z)X(z)$, which in this case is (I'm pretty sure)
$$
Y(z)= \frac{z}{z-0.4}\cdot 4 \frac{z}{z-0.25}\\
\begin{align}
\frac{Y(z)}{z} &= \frac{4z}{(z-0.4)(z-0.25)}=\frac{C}{z-0.4}+\frac{D}{z-0.25}\\
&=\frac{32}{3} \frac{1}{z-0.4} - \frac{20}{3}\frac{1}{z-0.25}\\
\end{align}\\
\therefore Y(z)=\frac{1}{3} \left[\frac{32z}{z-0.4} - \frac{20z}{z-0.25}\right]\\
\therefore y[n]= \left[ \frac{32}{3}(0.4)^n - \frac{20}{3}(0.25)^n \right] u[n]
$$
which is almost the same, but not quite.

Where did I go wrong? Also, in the second approach, how do you take into account the initial condition of $y[0]=0$? And thirdly, in the first approach, is the first term of $y[n]$ multiplied by $u[n]$ or not? If so, how do you show this?

(for anyone wondering, this is not homework for a course, just self study...)

As Matt L pointed out, your initial condition is wrong. y[0] is NOT 0, it's 4. This is obvious from the difference equation. Once you use y[0] = 4, both results come out to be the same.


Average Power Spectral Density of PAM signals I am reading through the PAM transmission scheme and about the power spectral density of the signals. Given that the Average Power Spectral Density of PAM Signals is:
$$
\Phi_{ss}(f)=\Phi_{aa}\left(e^{j2\pi ft}\right)\frac{\lvert G(f)\rvert^2}{T}
$$

(We consider that the amplitude coefficients of PAM $a[k]$ is a cyclo stationary discrete stochastic process and hence compute average of a PSD).

There are three remarks made after this:

It becomes obvious that the minimum bandwidth of the pulse $g(t)\overset{\mathcal F}\rightarrow G(f)$ has to be at least $1/T$.
The PSD $\Phi_{aa}\left(e^{j2\pi fT}\right)$ of the discrete time sequence $a[k]$ is periodic with period $1/T$.
If $G(f)$ has smaller bandwidth than $1/T$ a full period of the data spectrum is not contained in the PSD of the transmit signal.

I am not sure if I have fully understood the remarks. Here is my understanding on the above:

Since we have taken an Average PSD, the PSD repeats over the time
interval $T$ and the bandwidth of the pulse is $1/T$. (I am not sure how that became the minimum bandwidth).
The PSD of the discrete sequence $a[k]$ is periodic since it's a cyclo stationary process and vary over a time period $T$.
If $G(f)$ would have had less bandwidth than $1/T$, that is the pulse has a higher frequency than $1/T$, the data spectrum would have smeared into each other leading to ISI.

Please let me know if the above understanding is correct/valid. Please add/correct if anything is missing.

Question 1:

If we have a PAM sequence like this:

$a[k]=[0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0]$,

then the PAM signal (Waveform) must hold frequency components which cover these fast transitions from 0 to 1 to 0 to 1. The oscillation speed (frequency) of the transitions are basically given by the time $T$ between the PAM symbols in the waveform.

If $G(f)$ has lower bandwidth than $1/T$, these fast transitions are not possible. All the frequency components which could oscillate at that speed are filtered out by $G(f)$.

Here a plot of the above Waveform, its PAM sequence $a[k]$ and the single pulses (when $a[k]=1$) which summed up create the waveform.

And another plot with the above Waveform and its highest frequency component. If the dashed cosine would oscillate any slower, it could not represent the Waveform anymore.

Question 2:

The minimum bandwidth of $G(f)$ is given by $1/T$, but it can be higher and almost always is (see root raised cosine filter and roll-off factor, which let's you control the bandwidth of $G(f)$). BUT, the PSD of a PAM sequence is only defined within $1/T$, and due to the Nyquist–Shannon sampling theorem it is repeated thereafter (aliasing).

Why is this important?

Look at the following plot, it shows the PSD of a PAM sequence (limited by $1/T$), the PSD of a PAM signal and the filter $|G(f)|^2$.

$G(f)$ filters the PAM sequence and the result is the PAM signal, but what happens outside $1/T$? $G(f)$ filters the aliased part of the PAM sequence.

Note: Not sure if the y-axis of the PSD is correct, but we are only interested in the bandwidth.

Question 3:

You are correct!

But there are some tricks (which e.g. are exploited by the (root) raised cosine filter), see my first plot and the single pulses. They obviously overlap in time, but at each symbol a[k] only one pulse is responsibly for all the amplitude of the waveform. All the other pulses related to other symbols are zero. See the Nyquist ISI criterion.


Detecting & Removing steps in signal Given the signal shown below, what is the best way to remove the steps and local maximas it contains. The signal contains some steps which can last up to 100 Samples before they return to about the same value as before the step (Marked with red circles). There are also some peaks which last only for one sample which should also be removed.

First Dataset: On this set, the median Filter works as expected (not shown here).

On the following Dataset, the Medianfilter with the same settings fails. To many samples are replaced by the median value.

What I've tried already:

Work with a moving average. Advantage: Small peaks have less influence and are nicely removed. Disadvantage: Long lasting steps remain in the signal.
Get the first difference of the signal (diff(X)) and correct the value in X if a big positive change is followed by a negative change immediately. Advantages and disadvatages are the same as above.

What are some good methods to remove/correct these unwanted values?

The problem with the moving average is that the average is not robust to the presence of the outliers - so you would need a very large window size to "dilute" the outliers.

Try a non-linear filter instead, like a median filter: Apply a median filter on your signal - you would need a window size of at least 300 samples. Compute the difference between the original signal and the median-filtered version. If the difference is above a threshold, replace the signal by the median-filtered version.

Here is some scilab code that attempts to implement this suggestion.  The results are plotted here; it seems to work nicely on the fabricated data.

function sm = smooth(x,len)
    sm = filter(ones(1,len)/len,1,x);
endfunction

N = 6000;
x = 0.1*rand(1,N,'normal');
y = cumsum(x);
y = smooth(y,100);

clf
subplot(211);
plot(y)

Njumps = 20;

jump_indices = round(rand(1,Njumps)*N);

jump_length = 0;

y2 = y;

for idx = jump_indices,
    y2(min(N,idx:(idx+jump_length))) = y(min(N,idx:(idx+jump_length))) + 1;
    //plot(min(N,idx:(idx+jump_length)),y2(min(N,idx:(idx+jump_length))),'r')
    jump_length = jump_length + 1;
end

plot(y2,'g');

filter_length = 10;
y3 = y;
for k = 1:N,
    y3(k) = median(y3(max(1,min(N,(k-filter_length/2):(k+filter_length/2)))));
end

plot(y3,'k')

subplot(212);
plot(y-y2);
plot(y-y3,'r');


IIR filter design in digital domain using the magnitude squared Does anyone have any good references for deriving parameters of an IIR Low pass/High Pass filter directly in the digital domain using the magnitude squared at the corner frequency?

I have been able to derive the parameters of a first order Low/High pass filter with $3\textrm{ dB}$ attenuation at the corner frequency i.e. calculating $k$ and $\alpha$ in:

$$H(z) = k\frac{\left(1+z^{-1}\right)}{\left(1-\alpha z^{-1}\right)}$$

My issue is that I distinctly remember deriving the parameters using a $6\textrm{ dB}$ attenuation at the corner frequency in a DSP course I have done previously but I have forgotten the trigonometric identiftes used to finish the derivation.

The general procedure is as follows:

Let $\omega = 0/\pi$ to calculate the gain term $k$ such that there is a $0\textrm{ dB}$ gain at $0/\pi$
Calculate the magnitude squared at the corner frequency to obtain a value for $\alpha$ in terms of the corner frequency.

The problem may be that it should be a second order filter or I am recalling the method for a band pass/stop filter but I'm not sure and it appears this method is not used very often except in the case of band pass/stop filters for parametric EQ.

I hope the question is clear and I will try to improve the structure with the responses so it will be useful for others. Any help will be appreciated.

To solve the case that you mentioned...

You have 2 variables to determine, so you need two relationships to resolve the two variables.  I'm going to use $k$ and $a$ as the variables to make this easy to type up.

$$H(z) = k\frac{1 + z^-1}{1-az^-1}$$

Start by considering the passband gain.  Use $f = 0$ for this.  Assume you want unity gain at $f_0=0$.

Assume:  $H(f_0) = 1, f_0 = 0$

Substitute $e^{i2\pi f/f_s}$ for $z$, $f_s$ is your sampling rate, set $f = 0$ and solve for $k$ to satisfy $H\left(f_0\right) = 1$.

From this you get $k = \frac{(1-a)}{2}$

Now work on the gain squared at your desired corner frequency ($f_c$) to determine $a$.

$H(f_c) = -3\textrm{ dB}$ (magnitude squared will be $-6\textrm{ dB}$ as you've stated)

We'll work with the magnitude squared at $f_c$ and set the gain to $1/2$ ($-6\textrm{ dB}$).

$\lvert H\left(f_c\right)\rvert^2 = \frac{1}{2}$

This time substitute $e^{i2\pi fc/f_s}$ for $z$.

To simplify the arithmatic you can solve this equation:

$$
\left(\frac{\lvert H(f_0)\rvert}{\lvert H(f_c)\rvert}\right)^2 = 2
$$

This eliminates the factor $k$.

You will end up with a quadradic relationship in $a$.  Solving for $a$ yields:

$$
a = \frac {1 - \sqrt{1-\cos^2\left(2\pi\frac{f_c}{f_s}\right)}} {\cos\left(2\pi\frac{f_c}{f_s}\right) }= \frac {1 - \sin\left(2\pi\frac{f_c}{f_s}\right)} {\cos\left(2\pi\frac{f_c}{f_s}\right) }
$$


Drift in PID controller I am using MATLAB Simulink's toolbox for Parrot Minidrones, which have an ultrasound an optical flow sensor built-in. Whenever I am trying to let the drone hover in position, it drifts into one direction (always the same direction). The strange thing is that the drone knows it is drifting, but nothing is happening to correct it:

The drift only seems to occur horizontally, vertically the drone is not drifting. An overview of the relevant part of the PID implementation can be found here. I have tried tweaking some values, but it does not seem to improve. Does anyone have any tips for what the issue might be?
Edit: graph of roll, pitch & yaw

Your controller controls pitch and roll. However, in the graph you showed us we only see x,y and z.
Could you add the plot for pitch and roll? I'm not an expert in aircraft control, but I think you can have a constant pitch and roll and still move in the x,y,z plane.
I'll edit my answer as soon as you give more info.
Edit : You roll seems stable, you could try adding an integral component to your pitch roll controller.


Block code error control with $d \ge c+u+1$ This is about error control coding formula

$$ d \ge  c+u+1$$

of a block code, with

$c$ being the number of correctable errors,
$u$ being the number of uncorrectable errors, and
$d$ being the distance of the code.

Per my understanding, this is applicable for all the block codes - linear or not. Now let us dissect the repeat code with respect to this formula.

Repeat code with $k=1$, $n=7$: Let us say we repeat each bit 7 times. $\mathrm 0$ is sent as $\mathrm 0000000$ and $\mathrm 1$ is sent $\mathrm 1111111$.  $d = 7$ in this case.  We conclude that all these combinations below are possible:

$$\begin{array}{|c|c|}\hline
 c & u \\ \hline
6&   0\\
5&   1\\
4&   2\\
3&   3\\
2&   4\\
1&   5\\
0&   6\\\hline
\end{array}$$

Now let us consider $c=3$ and $u =3$.  With respect to the repeat code we are discussing above, what does it really mean?
Can someone kindly explain with examples ?
How does the explanation change with one or two more of the above combinations?

You have to realize that $d$ is specifically the minimum distance in your code, i.e. the distance between the to code words that are closest together.

And $u$ is not exactly the number of "uncorrectable errors", but the maximum number of errors that can definitely be detected, but not corrected correctly.

Note that "distance" here is the hamming distance, i.e. the number of bits you'd have to flip until you end up with another codeword.

  We conclude that all these combinations below are possible: ...

False. The first three lines ($c\in\{6,5,4\}$) are impossible.

Also, your table doesn't make much sense. Given a code, you'd always correct as many errors as possible¹, so for a given code, $c$ is given, and $u$ is given. The fact that

$$\begin{align*}
c &=& \#(\text{detectable bit flips, which ARE correctable})\\
  &=& \#(\text{total detectable bit flips}) &- \#(\text{detectable bit flips, which AREN'T correctable})\\
  &=& (d-1) &- u\tag{*}
\end{align*}$$

is just your definition of $c$.

A $k=7$ repetition code can never correct 6, 5 or even 4 bit errors; I find "you can't flip more than $\left\lfloor \frac d2 \right\rfloor$ bits until it looks more like another codeword than like the original" highly intuitive, so the only way I can explain this is by example:

take your 0 code word, and flip 4 bits:

0 0 0 0 0 0 0
0 1 1 1 0 1 0

Now, that looks much more like the 1 keyword, doesn't it? In fact, 0 1 1 1 0 1 0 has a hamming distance of 4 to the 0 code word, and a hamming distance of 3 to the 1 code word, so it's closer to 1, and will correctly be identified as error, but falsely get corrected to 1.

Looking at this, any error is detectable, as long less than the minimum distance bits are flipped (for any code, pretty logical, isn't it?). So,

$$c_\text{max} = d -1 \text.$$

Now, your formula simply "happens" when you put in your definition of $c$ and $u$, see $(*)$ above.

There's not much to explain here: for repitition codes $d=k$, and you can always correct $\left\lfloor \frac d2 \right\rfloor=\left\lfloor \frac k2 \right\rfloor$ bit errors and detect $d-1=k-1$ errors.

One can't explain your other combinations, as they are basically wrong for any maximum likelihood decoder. You can of course work "worse" than possible, but especially for the  repitition codes, summing up the received bits and deciding whether the result is larger than $\frac k2$ gives you $c=\left\lfloor \frac k2 \right\rfloor$ and thus, $u=(d-1)-c=(k-1)-\left\lfloor \frac k2 \right\rfloor =\left\lceil \frac{k-1}{2}\right\rceil $, and, in your $k=7$ case, $c=u=3$.

¹: point is that there's many codes where you can't, with reasonable effort, correct as much correctable errors as the code would allow you too, that means that $c \le \left\lfloor \frac d2 \right\rfloor$ in general, and that's why $d\ge c+u +1$ and not $d = c+u +1$. However, repetition codes are so trivial, we can truly Maximum-Likelihood-decode them and correct $\left\lfloor \frac d2 \right\rfloor$ errors.


Are all LTI systems invertible? If not, what is a good counterexample? I have been trying to figure this out for a while now. Everywhere I have looked I could easily find examples of invertible LTI systems, but I could not find any counterexamples. Can anybody shed some light on this for me?

You need to define what you mean by "invertible". Do you mean invertible by a causal and stable system? If yes, then any system that is not minimum-phase is not invertible (because the inverse system can't be causal and stable).

Example of a system that cannot be inverted by a causal and stable system: a simple delay $y(t)=x(t-T)$, $T>0$, could only be inverted by a non-causal system.

If by "invertible" you mean a system that can be inverted by a system that is possibly not causal and/or not stable, then it's still straightforward to find systems that can't be inverted according to this broader criterion: just take a system with one or more zeros in its frequency response. The information about the input signal at the frequencies where the system's frequency response is zero cannot be recovered by any other system.

Example of a system that cannot be inverted by any system: an ideal low-pass filter cannot be inverted by any system, because any information above the filter's cut-off frequency is lost and cannot be recovered.

EDIT: the latter definition of invertibility (where causality and stability of the inverse system are ignored) is just the definition of an injective mapping. In terms of systems this means that no two distinct input signals generate the same output signal. This category is what is also described in Laurent Duval's answer. Note, however, that in practice we usually prefer inverse systems that can also be realized, i.e., we require stability and causality, and then the first definition of invertibility given above is appropriate.


Fitting discrete curve to cosinus I have a curve which the asymptotic behavior is the one of a cosinus. I look at this curve on all integers.

I know that the asymptotic curve has this behavior : $$ cos(a*x+b) $$

My question is : how can I determine $a$ and $b$ ?

Indeed I only have samples of the curve, and the samples that I have are not exactly of the form of the cosinus as the cosinus is the asymptotic behavior of my curve.

I did the following : I took samples between $[N; N+P]$ with $N$ large enough to be sure that I'm ""in"" the asymptotic zone.

I did a FFT on this zone and noted the frequencies where a Dirac seemed to appear. I . (I chose P large enough to have enough sample to have a good frequency resolution).

Then I made a function which fft is just 2 pure dirac.

But when I plot my curve and the supposed asymptote, I have this (in blue the initial curve that I evaluate only on integers and in red the function I made with the two diracs).

We can see that we have a slow drift between the two curves (I'm almost sure that if I look to samples far away the two curves will be in opposite of phase at a moment).

My question is : is it possible to determine an asymptotic behavior of samples in my case ? And how to do it ?

As pointed out by Dan Boschen, your signal appears to have a frequency, which does not lie exactly on an FFT bin. Given that the signal contains only a single cosine, you can get a better estimate of its frequency by zero-padding the signal. Have a look at the code below, which reproduces your problem:

Fs = 450
f = 75.5

T = 1

t = np.arange(0, T, 1/Fs)

x = np.cos(2*np.pi*f*t)

plt.figure(figsize=(10,10))

plt.subplot(221)
plt.plot(t, x)
plt.xlim((0, 0.1))

f1 = np.linspace(0, Fs, len(x), endpoint=False)
plt.subplot(222)
X1 = abs(np.fft.fft(x))
plt.plot(f1, X1, '-x')
plt.xlim((60, 90))
plt.axvline(f)

f2 = np.linspace(0, Fs, 32*len(x), endpoint=False)
plt.subplot(223)
X2 = abs(np.fft.fft(x, len(f2)))
plt.plot(f2, X2, '-x')
plt.xlim((60,90))
plt.axvline(f)

plt.subplot(224)
plt.plot(f2, X2, '-x')
plt.xlim((72,79))
plt.axvline(f)

iMax1 = np.argmax(X1)
iMax2 = np.argmax(X2)
print ("Estimated frequency of the Cosine without Zero-padding: f=%f" % (Fs-f1[iMax1]))
print ("Estimated frequency of the Cosine with Zero-padding: f=%f" % f2[iMax2])

The top-left figure shows a part of the sine wave. The top right shows its fft (zoomed into the interesting frequency). As you can see, the FFT shows a significant values at two different frequency bins. This indicates that the actual frequency is between both bins. To find the correct frequency, we can apply zero-padding. Bottom left shows the zero-padded FFT, showing more samples in the frequency domain. The bottom-right shows a zoomed-in version and we see, that the peak of the FFT is exactly at the correct frequency f=75.5Hz. For more information about zero-padding, you can refer to an article I wrote about when you can apply zero-padding, when it wont help etc.

The textual output of the code above is

Estimated frequency of the Cosine without Zero-padding: f=76.000000
Estimated frequency of the Cosine with Zero-padding: f=75.500000

As we see, if we take the maximum of the frequency bin of the non-zeropadded FFT the frequency is estimated to be at 76Hz. However, if we consider the zeropadded FFT, we get the correct frequency estimate.


what is the differnence between 8x8 dct coefficients and 25 x 25 dct coefficients? As I knkw the DCT is discrete cosine transform that mostly used in time domain analysis to frequency domain. Also 8x8 is used in nomaly in img processing. I don't know why that is defined by 8x8 .anyway 8x8 is used in processing.

BTW, I am want to know what is the difference between 25x25 and 8x8 in dct?
Is this just for quality?

There exist different variants of discrete cosine transforms. They mostly vary on a phase shift. The most commonly-used type is known as type-II. Since the works of K. R. Rao et al. in the seventies, it was shown to possess near-optimal properties.

The DCTs are invertible, possess certain symmetry properties that make them useful and practical in sound (signal) or image processing. They can technically be applied to any patch size 13, 8x8, 25x32... Power-of-two patches benefit from faster and more optimized computations. 8x8 are standard for natural image due to:

a size small enough so that image chunks are stationnary enoguh to be quite-well represented by cosine waves
a small memory footprint for digital implementation

Now, it is more common to find 4, 8, 16, and 32 sizes in image and video coding, and even non-square blocks (4x8) with close to DCT transforms.

The quality is not a real question for the DCT, since it is invertible. It is more related to a match between your data, the transform on the processing (e.g. compression). You might have data and applications where a 25x25 DCT would work better. Then you will probably hit some difficulties for optimized implementations.


Counter design for audio signals I am designing a schematic in which a Counter is to be increased whenever a particular band of frequencies is received at the Audio input.
We are receiving the audio input as electrical Signal and I have already designed the bandpass filter with the required passband. How can I design the Counter for the same purpose?

def doppler(x):
    """
    Parameters
    ----------
    x : array-like
        Domain of x is in (0,1]

    """
    if not np.all((x >= 0) & (x <= 1)):
        raise ValueError("Domain of doppler is x in (0,1]")
    return np.sqrt(x*(1-x))*np.sin((2.1*np.pi)/(x+.05))

def blocks(x):
    """
    Piecewise constant function with jumps at t.

    Constant scaler is not present in Donoho and Johnstone.
    """
    K = lambda x : (1 + np.sign(x))/2.
    t = np.array([[.1, .13, .15, .23, .25, .4, .44, .65, .76, .78, .81]]).T
    h = np.array([[4, -5, 3, -4, 5, -4.2, 2.1, 4.3, -3.1, 2.1, -4.2]]).T
    return 3.655606 * np.sum(h*K(x-t), axis=0)

def bumps(x):
    """
    A sum of bumps with locations t at the same places as jumps in blocks.
    The heights h and widths s vary and the individual bumps are of the
    form K(t) = 1/(1+|x|)**4
    """
    K = lambda x : (1. + np.abs(x)) ** -4.
    t = np.array([[.1, .13, .15, .23, .25, .4, .44, .65, .76, .78, .81]]).T
    h = np.array([[4, 5, 3, 4, 5, 4.2, 2.1, 4.3, 3.1, 2.1, 4.2]]).T
    w = np.array([[.005, .005, .006, .01, .01, .03, .01, .01, .005, .008, .005]]).T
    return np.sum(h*K((x-t)/w), axis=0)

def heavisine(x):
    """
    Sinusoid of period 1 with two jumps at t = .3 and .72
    """
    return 4 * np.sin(4*np.pi*x) - np.sign(x - .3) - np.sign(.72 - x)
In [5]:
x = np.linspace(0,1,2**11)
dop = doppler(x)
blk = blocks(x)
bmp = bumps(x)
hsin = heavisine(x)
In [6]:
fig, axes = plt.subplots(2, 2, figsize=(10,10))
ax1 = axes[0,0]
ax2 = axes[0,1]
ax3 = axes[1,0]
ax4 = axes[1,1]

ax1.plot(x,dop)
ax1.set_title("Doppler")

ax2.plot(x,blk)
ax2.set_title("Blocks")

ax3.plot(x,bmp)
ax3.set_title("Bumps")

ax4.set_title("HeaviSine")
ax4.plot(x,hsin)

for ax in fig.axes:
    ax.tick_params(labelbottom=False, labelleft=False, bottom=False,
                   top=False, left=False, right=False)

fig.tight_layout();


Calculating a spectral envelope I'm trying to calculate the envelope of an arbitrary spectrum but I'm slightly lost as to the best way to do this.

I've looked at cepstral enveloping but its very fiddly and I can't seem to get it to give the right results.  I've also tried LPC but again it eems to miss a lot of the minor peaks while enveloping the major peaks perfectly.

I'm beginning to wonder of the best way would be just to run a sliding window across each frequency bin in an FFT returned magnitude spectrum and select the maximum value within that window.  This way I'd simply pull out all the peaks (at whatever granularity I'm after) and then I can simply plot a curve through the resulting points.

Is there a correct way to calculate the spectral envelope?

Cepstral windowing is the most obvious way to extract a spectral envelope, you have a description here.

  I've looked at cepstral enveloping but its very fiddly and I can't
  seem to get it to give the right results.

In the linked description, try lowering the lowpass "cut-off" sample (nc) until you retain the amount of detail you want. You can get useful implementation examples here.

  I've also tried LPC but again it eems to miss a lot of the minor peaks
  while enveloping the major peaks perfectly.

Isn't that the point of extracting the spectral envelope? By LPC are you referring to the Linear Prediction Spectral Envelope?


I have the scatter plots for a Hamm and Tukey(0.25) windowed data. The complex data plotted from Hamm windowed forms a pattern, why? The Hamm windowed complex data forms some sort of a quantized pattern whereas in the same code if tukeywin(0.25) is applied normal noise is observed. Any idea why ?

Assuming we're looking at plots of the time domain I&Q, I'll guess that your multiplication operation is of limited bit width and is causing the quantization.

If you look at the Hamming Window in the time domain, only the central sample(s) are weighted at 1.0, with all the other weights down along a relatively smooth curve.
https://en.wikipedia.org/wiki/File:Window_function_and_frequency_response_-Hamming(alpha_%3D_0.53836).svg

If you look at the Tukey Window most of the samples are weighted at 1.0, with some samples on the edges of the window weighted along a curve. https://en.wikipedia.org/wiki/File:Window_function_and_frequency_response_-Tukey(alpha_%3D_0.5).svg


Fourier transform diagonalizes time-invariant convolution operators I got the following paragraph from the book "A wavelet tour of signal processing" chapter one, page 2.

The Fourier transform is everywhere in physics and mathematics because
it diagonalizes time-invariant convolution operators. It rules over linear time-invariant signal processing, the building blocks of
which are frequency filtering operators.

How is it (illustrated) formulated mathematically?

TL;DR: Fourier is the logarithm of the convolution.
Doing a convolution (linear time-invariant operator) is a quite involved processing in the original (primal) domain: it involves intricate sums and products of an operator $h$ and signals $s$ at different indices:
$$ \sum_k h[k]\times x[n-k]\,.$$
With a Fourier transform, you transform the whole information for the operator $h$ and data $x$ in the primal domain onto a novel dual (frequency) domain, in $\hat{h}$ and $\hat{x}$. Then the result of convolution becomes:
$$  \hat{h}[m] \times \hat{x}[m]$$
which is a simple index-by-index product. Hence the term "diagonalization": from products of two indices, $[k]$ and $[n-k]$, Fourier turns them into a product of a single index $[m]$ times  $[m]$.
To make the above a bit simpler, Fourier is to convolution what logarithms are to multiplication. Let us forget the carry for a moment. Doing a multiplication is complicated, especially for large numbers: one has to multiply thousands to tens and add them to products of hundreds.  Logarithms turn multiplication into adds: you (almost) have to add figures in the same position (lest the carry), instead of combining figures in different powers.
For ages, people have used logarithm tables or logarithmic sliding rulers as analog computers. Logarithmic tables and rulers have been very important in days before computers, to allow computations in engineering, astronomy, etc.
Fourier principles, tools and fast librairies are the nowadays analogs of logarithms, that helped shaped our digital world. The FFT algorithm is often considered one of the most important algorithms. Indeed, there exists FFT based multiplication of large numbers.


How to do image reconstruction from array data? i have 15 set of data which is the s-parameter, consist of 2 column (1001x2) each set. The first column is frequency and the second column is magnitude in dB. I have transform the magnitude from (dB) into value and frequency into time. But i have no idea how to transform these 15 sets of data into an image by using Matlab. Can anybody help me?

You will not be able to reconstruct the image without the phase information.  The frequency magnitudes alone are not sufficient.


Determine the haar wavelet function is linear or nonlinear I'm looking for analytical justification of linearity or non linearity of the wavelet transform with the real Haar mother wavelet function.

I have googling already. But I can't find and understand the meaning

Linearity can play at different levels in this context:

The Haar wavelet functions are piece-wise affine (constant), not linear (but sometimes people call affine "linear").
Any standard wavelet decomposition is linear, as it decomposes a signal into a sum of coefficients times a wavelet. So this is not specific of the Haar wavelet.
Wavelet approximations are often nonlinear, as they consist in choosing the $K$ best (largest) coefficients to represent the signal. And choosing the $K$ largest does not obey the linearity principle.

So I cannot answer directly, by lack of context.


In the Sampling Theorem, why are the image frequencies at n*fc not a problem? In this example I have been working through, we first look at the situation when fs > 2fc and then the situation when it isn't:

In the example, the frequency responses of a sampled sinusoid at fc = 5 kHz are shown. In the first instance, the sampling rate is fs = 16 kHz, and in the second it is fs = 8 kHz.

I understand the concept of aliasing and I can see how it has occurred in the second example - an alias has crossed over fc. What I don't quite get is why the situation in the first instance is acceptable - there are still aliases in the signal. Is it simply that we can window (low-pass filter) the frequency response to fs/2 as we know this is the highest frequency that can be resolved?

The "alias" frequencies you see appearing above the Nyquist limit $f_s/2$ do not exist in a proper interpretation of the digitised signal; they are appearing due to incorrect interpolation of the samples during reconstruction. If you correctly interpolate the signal from the samples, those frequences will not be there.

Remember, a key condition of the Nyquist–Shannon sampling theorem is that it can perfectly record the input signal only so long as it contains no frequencies above $f_s/2$. A corollary of this is that the output signal you produce from correctly interpreting the samples must not contain such frequences either. The core of the theorem is that if you sample a signal under this condition, and then plot a sequence of points corresponding to each sample, there will be only one curve you can interpolate that both passes through every point and also contains no frequencies above $f_s/2$; that curve is the same as the original input.

The difficulties and details of adhering to this condition and properly interpolating the output are what cause problems like the one you're seeing.

In the case of your output, you're probably reconstructing the signal using piecewise interpolation, producing a "stair-step" output that looks like this:

(Image credit: https://wiki.xiph.org/Videos/Digital_Show_and_Tell)

That's certainly one way of doing it, but it (or at least it used alone) is the wrong way. Such an interpolation clearly has a lot of frequency information above $f_s/2$ (the alias frequencies) and is significantly different the original input signal meeting the conditions above.

It is, however, often an intermediate signal during the interpolation process; it turns out that one way of doing the interpolation fairly accurately is to to generate a signal using piecewise interpolation and then filter out all components of it above $f_s/2$, which will leave us with the original waveform.

You can do this in your first diagram by simply erasing everything above $f_s/2$, leaving only the original 5 Khz frequency you sampled.

For a more in-depth explanation of this, I recommend Monty Montgomery's, "D/A and A/D | Digital Show and Tell" video (also on YouTube). There's also a text version if you're short on time, but the video demonstration is signifantly more clear. I would consider this almost mandatory viewing for anybody involved in any details of digital sampling, audio or otherwise.


Conversion from stationarity to non-stationarity Is there any way to convert a non-stationary signal to a stationary one, perform operations on it meant for a stationary signal and then convert it back to the non-stationary one?

You've asked, if it's possible to convert a non-stationary signal into a stationary one and vice versa; the answer is yes it's possible to some extend as outlined below.

Indeed it forms the basis of a family of transform based signal compression algorithms. Specifically, for transform based signal compression, a block of time domain samples $x=\left[x_0,x_1,...,x_{N-1} \right]$ are taken, and transformed into a block of frequency domain coefficients $c=\left[c_0,c_1,...,c_{N-1}\right]$ practically via a DCT (discrete cosine transform) or similar and are re-quantized for furher processing.

Note that the time domain signal in the block is assumed to be WSS and correlated, as each sample $x[n]$ has the same variance $\sigma_x^2$, whereas the frequency domain coefficients are non-stationary with each coefficient $c[k]$ having a different variance, $\sigma_{c_k}^2$. This change in variance results in change in bit-resolution of the quantizer, hence compression. Furthermore, the compression efficieny is a direct consequence of the strength of the correlation in the time domain samples.

This non-stationary frequency transform coefficients $c[k]$ are made use of opimal quantization of those coefficients with larger variances (typically the lower frequency coefficients) assisted by the human visual system perceptual characteristics as well. That's a lossy compression. However, you don't have to compress anything, if all you need to get is non-stationary signal.

So this provides an answer for your question, but whether a sutiable method for your financial signal exists is up to you.

Have a look at the following books:

A Primer for Financial Engineering: Financial Signal Processing and Electronic Trading_AKANSU

Financial Signal Processing and Machine Learning_AKANSU


filtering in frequency domain I am performing decorrelation of two channels of audio using all pass filtering. The two all pass filters for two channels are of very high order (>1000 but with not all coefficients non zero). I have to perform this filtering in frequency domain. I get input samples in 2048 frame size (48 K sampling) with 50 % overlap Hann windowed.

So to apply this filter (a small intermediate matlab snippet for one of the channels)

freq_resp_filter=freqz(filter_APF,4096,'whole');

fft_input_sig=fft(input_sig_frame,4096);

filtered_sig_fft=freq_resp_filter.*fft_input_sig;

filtered_sig_ifft=ifft(filtered_sig_fft); % done for testing only

final_2048_samples=filtered_sig_ifft(1:2048); % throw away last 2048 samples to
remove aliasing by fft convolution

and then i recombine this final_2048_samples with previous frames with overlapping.

When I listen to the output audio, I hear intermittent noise with signal,
the noise is absent when I apply the same filter in time domain.

my question is:
1. Is this implementation correct (I have to perform this in frequency domain! and I can not do overlap add/save as this block is a part of a system, and I have no control over frame size, so every operation, 2048 samples should be given as output for further processing. )

Regards,
Kd

Here is a visual answer why it is not working. The overlap is not even symmetrical in the last step so the scheme can't work.

Even if you make it symmetrical by making the filter zero-delay, the amount of stuff that gets summed together varies by time, so it can't be right.


How can I estimate the spatial resolution of consumer grade electronics? I have been looking at stereo vision applications, and one working example involves using two USB web cams, mounted together and calibrated.

Is there a way to estimate the theoretical spatial resolution of such a set up?  The goal is to avoid purchasing and testing blindly but rather to be relatively confident ahead of time.

I do understand that the EXACT resolution will probably require a test setup due to lensing, etc.

What information would be needed to make this determination?  Would the minimum focal length and the maximum resolution be sufficient?  Pixel density obviously plays some role in this, and is advertised on most digital cameras (which have a much higher maximum resolution), if I was attempting to compute an estimate, how should I incorporate this?

As an example consider the following setup:

a 4" minimum focal length (MS HD Webcam focal length)
a maximum resolution of 1920x1080
object is on the order of a couple of centimeters
details on the order of tenths of millimeters
the distance between cameras and objects should be smallish, I was considering something like 15 cm to 25 cm, flexible based on the answer here
the distance between cameras, the number of cameras can be flexible but should stay in the smallish range ... if I need $10,000 in cameras, it's a bad solution

Consider the following figure:

The cameras are separated by baseline 2b - say 10cm, so b=5cm.  Your focal length of 4" sounds looong (about 100mm) for a webcam, but hey, let's go with it.  f=100mm (or also 10cm!).

From $xr-xl = \frac {2bf}{z}$ an object at 10m away will have a pixel difference between the two cameras ($xr-xl$) of 0.1*0.1/10 = 0.001m.  So, if your pixels are 5 micron square $xr-xl$ is about 200 pixels.

Which is a lot - large disparities are expensive to search for.

With a more normal focal length of say 5mm, your disparity is more like 10 pixels.  11m is about 9 pixels, so your distance resolution is about 1m @ 10m.

If you have a specific camera in mind you can run the numbers for your particular geometry.


What is the name of this filter? I am having this simple function which acts on two state variable $s_t$ (states being $\{0,1\}$, boolean in sense), and I am evaluating it in discrete time steps. The function is $f_t=1$ if $s_t$ has been $1$ for last $W$ samples (some fixed window size), and $f_t=0$ otherwise. So basically this will set the output to active if the input signal has been active for certain amount of time/samples. Is there some name or category for this? I thought it is latching but according to definitions on web it seems latching is more about delaying the signal, instead of asking "was the signal delayed on input".

There's no known name for this kind of filter in the signal processing world, as far as I know – and what good would such a name be, anyway?

So, I'd describe it in a signal processing context as "discrete state moving window sum with threshold", since it's a non-linear filter (and not much of a signal filter, if you ask me!). Simply because that's how I'd implement it: Sum up the last $W$ states, and output a 1 iff the result is $W$.

If you leave the realm of signal processing and think more of the basic electronic circuits that implement that, or by the purpose: You could call it a watchdog, or a digital anti-ringing, or a self-resetting cycle counter with a threshold, for example.

Your approach with the latch is from the latter category of names; "latch" is not anything that I'd call a filter in a signal processing context. Especially, the functionality you describe is not a latch (alone).


What are the tradeoffs for more coefficients for fractional delay filter? Is there a rule of thumb for how many coefficients to choose for a fractional delay filter?  I figure there might be a trade off based on the signal frequency as well as my sample rate (this is running on an FPGA), but I can't seem to find a good rule of thumb.

More taps = more resources (multipliers and registers)
More taps = longer transient response, in some applications it can be a dealbreaker. It could cause significant pre-echo which can be problematic in some applications such as ultrasound time-of-flight diffraction.
More taps = better frequency response, your fractional delay will have a unity gain for a wider frequency range. Same for group delay.


Visualisation of bandwith of filter with sampling I'm having some trouble visualising (and understanding) the filtering of a sampled signal.

If I understand correctly, after sampling a signal, it is desired to have it go through a lowpass to retrieve the original signal again. The literature states the following below.

$\omega_0 < B < \omega_s - \omega_0$

Where $\omega_s$ is the sampling frequency and B is, what I assume, the bandwith of the filter.

I do not understand what is meant by the above expression and the literature is kind of crap at explaining it.

I was hoping someone can visualise this for me.

$B$ is the cut-off frequency of the anti-alias filter needed before sampling; $\omega_0$ is the signal bandwidth, and $\omega_s$ is the sampling frequency.

$B$ can also be the cut-off frequency of the (ideal) reconstruction filter (band-limited interpolator).

If this is the case, this suggests that the reconstruction filter is an ideal low-pass filter with cut-off frequency between the signal bandwidth and the next overlapped frequency component ($\omega_s - \omega_0$). The sampled signal contains repetitions of the original specturm every $\omega_s$, so when reconstructing, you want to get rid of all these repetitions and keep only the spectrum centered at 0.

I don't know what literature you are referring to, and why is it so crappy.


Determine whether the system is a FIR or IIR by looking the transfer function I have the following system:
$$
y[n]=\frac{1}{3}(x[n+1]+x[n]+x[n-1])
$$
After the Z-Transform we get

$$
\frac{y[z]}{x[z]}=\frac{z^2+z+1}{3z}
$$
which is of course the transfer function of the system. Now, I know that it would be a FIR if $\sum h(z) < \infty$ with $h(z)$ the impulse response. But As far as I know, it should be clear if it is FIR or IIR by just looking the transfer function. How can I know that?

The forward $z$-transform is written with negative powers on $z$. So, $z$ transform of your difference equation would be:

$$
Y(z) = \frac{1}{3}\left(zX(z) + X(z) + z^{-1}X(z)\right)
$$

yielding a transfer function of:

$$
H(z) = \frac{Y(z)}{X(z)} = \frac{1}{3}(z + 1 + z^{-1})
$$

The easiest way to determine whether a filter is IIR or FIR is to identify its pole locations. For FIR filters, there is a rule for this that is based on the structure of the impulse response:

If the system is causal (i.e. it is zero for all $n < 0$), then it is FIR if all of its poles are located at the origin ($z=0$).
If the system is anticausal (i.e. it is zero for all $n \ge 0$), then it is FIR if all of its poles are located at $z=\infty$.

Your system is neither causal nor anticausal. Instead, it has a causal component ($\frac{1}{3}(x[n] + x[n-1])$) and an anticausal component ($\frac{1}{3}x[n+1]$). So, it can be helpful to look at each component separately:

$\frac{1}{3}(x[n] + x[n-1])$: The $z$-transform of the causal portion of the transfer function is $\frac{1}{3}(1 + z^{-1})$. This expression $\to \infty$ as $z \to 0$, so it has a pole at $z = 0$. This meets the criterion for an FIR response.
$\frac{1}{3}x[n+1]$: The $z$-transform of the anticausal portion of the transfer function is $z$. This expression $\to \infty$ as $z \to \infty$, so it has a pole at $z = \infty$. This meets the criterion for an FIR response.

This would imply that you gave a difference equation for a non-causal FIR filter.


Homodyne concept of receiver I would like to understand a concept of homodyne architecture of I-Q demodulator and how the wave looks like at each step.

It is a bandpass Signal:

Step 1: multiplication with cos and sin wave.  I represent cos and sin wave like Dirac function. For cos is even, sin is odd.

Step 2: filter with LPF

Step 3: summation

Will be it correct?

You can also look at it in the time domain. Your received signal is

$$s(t)=s_I(t)\cos(\omega_c t)+s_Q(t)\sin(\omega_c t)\tag{1}$$

where $s_I(t)$ and $s_Q(t)$ are the in-phase and quadrature components, and $\omega_c$ is the carrier frequency (in rad). We have to assume that at the receiver the carrier frequency and phase are known (coherent demodulation), or can be retrieved somehow.

After demodulation we obtain two signals:

$$\begin{align}s(t)\cos(\omega_ct)&=s_I(t)\cos^2(\omega_ct)+s_Q(t)\sin(\omega_ct)\cos(\omega_ct)\\&=\frac12\big[s_I(t)(1+\cos(2\omega_ct))+s_Q(t)\sin(2\omega_ct)\big]\tag{2}\end{align}$$

and

$$\begin{align}s(t)\sin(\omega_ct)&=s_I(t)\cos(\omega_ct)\sin(\omega_ct)+s_Q(t)\sin^2(\omega_ct)\\&=\frac12\big[s_I(t)\sin(2\omega_ct)+s_Q(t)(1-\cos(2\omega_ct))\big]\tag{3}\end{align}$$

The low pass filters filter out all components at $2\omega_c$, so at the outputs of the filters we're left with

$$r_I(t)=\frac12s_I(t)\tag{4}$$

and

$$r_Q(t)=\frac12s_Q(t)\tag{5}$$


Applying zero-phase filtering in the frequency domain - also works in real-time? In my application a filtering in the frequency space is desired. I designed two high-pass filters: $h_1$  Butterworth IIR filter, and $h_2$ FIR. I designed them using Python's SciPy functions.
Using signal.freqz I obtained the frequency response $H_1(k)$ and $H_2(k)$.
I now want to filter signal in a real-time, buffer-wise. I read buffers of certain size, say $N=2048$, and do a DFT and obtain $F(k)$. I then want to multiply the Fourier coefficient with the two filters frequency response: $G_i(k) = F(K) H_i(k)$ for $i=1,2$. Then I do inverse DFT to get the filtered signal.
(Remark: in practice I use overlap/add in the buffering, and so recovering back the filtered signal in time domain requires windowing and is more complicated.)
Now for the questions:
Multiplying $F(k)$ by $|H_1(k)|$ seems to work, even in real-time. But multiplying by the absolute value is like performing zero-phase filtering, which by theory is non-causal. However, since I apply it buffer-wise, this seems to be the loophole. Am I correct?
I also filtered using $F(k) H_2(K)$ and in Python it works just fine. However, giving the list of coefficient $\{ H_2(k) \in \mathbb{C} \mid k=0,...,N/2 \}$ for a colleague, for her to filter in another application/software, the filtering didn't seem to work. The filtered signal was a total mess (the unfiltered signal is a short speech by a woman). I verified with her that she applied the filter with the "mirror conjugate" completion so $\tilde{H}_2 = (H_2(k), \overline{H_2(N/2+1-k)}) \in \mathbb{C}^N$. Is there a theoretical reason for it not to work? (so we can narrow the reason to a bug in her code)

It's not a loophole.  A causal linear phase filter is identical to a non-causal zero phase filter in series with a delay of half the filter length.  If you check the results of your filter with an impulse input, you will probably see such a delay.
So it's not really real-time if you count the delay, but is often close enough to seem real-time if the filter's impulse response and/or block size is short relative to your real-time requirements.


DCT and mean difference of an image If we subtract the mean of image from the DCT coefficients of an image

what does it signifies?
and what type of distribution does it follow?

Example in MATLAB:

i=imread('cameraman.tif');
c=mean2(i);
d=dct2(i);
f=c-d;

what does f signifies?

You can investigate that if the size of the DCT block (which is the whole image in here) is $x\times y$, then

c=d(1,1)/sqrt(x*y);

That is, the DC element of the DCT matrix (i.e. d(1,1)) is actually the mean of the image scaled by $\sqrt{xy}$.

So you are basically subtracting a scaled version of the DC component form all coefficients, which makes not much sense at the moment, except if you provide more context.


Confusion regarding wav to MFCC conversion I have to compare two .wav files.

I am converting the wav file to an MFCC file. Why am I getting different MFCC matrices when using different Matlab toolboxes (one of them is giving a matrix having 20 rows while the other two are giving matrices having 13 rows, and even those two matrices are not same)?

Why are the two matrices different? I have heard that there are 13 MFC coefficients. Does it have something to do with the number of rows?

There are many variations in the ways MFCCs can be extracted. The most common ones are:

STFT parameters like frame rate, FFT size, FFT window.
Filter bank parameters like number of bands, whether filter normalization is used, and whether the filters are applied to the magnitude or power spectra.
Subtle variants in the definition of the discrete cosine transform.

And it is actually not correct to say that there are 13 MFCC - it is indeed very common to keep the 13 first coefficients, but you could very well use a larger cosine transform matrix and get more coefficients. The farther you go, the less information the coefficients will capture (since they will correspond to dimensions of decreasing variance - assuming the discrete cosine transform is a good approximation of a PCA on the mel coefficients data, which it is).

Check this to learn more about all the subtle differences between several "famous" MFCC extraction toolboxes.

In the end, such or such combination of MFCC extraction parameters might give you a difference of at most one or two percentage points at whatever task you are using the MFCC for (for example speech recognition), but there are more interesting things to try to tweak and improve than these parameters.


Design Filter with Arbitrary Magnitude and Phase Response in Matlab Does anyone know of a way to design a filter with arbitrary magnitude and phase response in Matlab?

for FIR, the simplest way is to inverse FFT you two-sided complex frequency response:

$$ h[n] = \operatorname{iFFT} \left\{ H[k] \right\} $$

(use a Kaiser window to window the FIR to a length that is acceptable to implement.)

i presume you have access to firpm( ) and/or firls( ).

so first divide your complex frequency response into real and imaginary parts:

$$ H_r[k] = |H[k]| \cos\left(\arg\{H[k]\} \right) $$
$$ H_i[k] = |H[k]| \sin\left(\arg\{H[k]\} \right) $$

send both $H_r[k]$ and $H_i[k]$ to the Parks-McClellan or Least Squares optimal FIR designer, but set ftype to 'hilbert' for $H_i[k]$.

sum the two resulting real-valued FIRs to get your result.

lastly, if you want to design an IIR to an arbitrary magnitude and phase spec, then that's a much bigger problem and consideration.  maybe use Prony to design to your impulse response $h[n]$ or use Matt's dissertation that i first learned of now, or use Greg Berchin's FDLS.  i dunno.  no quick solution for you there.


Understanding Analytic Signals I know that a function $f_a\in \text{L}^2(\mathrm{R})$ is an analytic signal if its Fourier transform is 0 for negative frequencies, i.e.

$$ \hat{f_a}(\omega)=0$ if $\omega<0 $$

We can characterize an analytic function by its real part $f = \text{Re}[f_a]$, and the Fourier transform of this function is

$$ \hat{f_a}(\omega)=\left\{\begin{matrix}
2\hat{f}(\omega) & \text{if }\omega\geq0 \\
0 & \text{if }\omega<0
\end{matrix}\right. $$

Now if we have a discrete signal $f[n]$ of size $N$, we can also compute $f_a[n]$ by setting the negative frequency components of its discrete Fourier transform to 0.

Now why is it that if we want $\text{Re}[f_a]=f$, then we must have:

$$ \hat{f_a}[k]=\left\{\begin{matrix}
\hat{f}[k] && \text{if }k=0,N/2\\
2\hat{f}[k] && \text{if }0<k<N/2\\
0 && \text{if } N/2<k<N
\end{matrix}\right. $$

Could that be explained?

It's not clear what you're asking, but I think you want to know why the DFT of a signal has extra conditions for it to real.

For the continuous case: the 2 is there because we need a function to be symmetric to have a real Fourier transform.

In the discrete case, the DFT obey the symmtery:

$$
X_{N-k} = X_{-k} = X_k^*
$$
i.e. conjugate symmetry (half the values are just complex conjugates of the other half). Thus $X_0$ and $X_{N/2}$ are real numbers whilst the remainder of the DFT is specified by $N/2-1$ complex numbers (hence the 2 in the second line).

Also, the answers to this question.

offer counterexamples to your first statement - that analytic functions have vanishing negative frequencies.


Simulation and theory of first-order low pass filter I am trying to teach myself some signal processing and control theory. I'm currently working through THIS tutorial. I've started an ipython notebook to keep my notes in.

Unfortunately, something is wrong with my very first simulation of a simple first-order low pass filter. The differential equation is
$$\dot{x} = \omega_0(u-x)$$
where $u(t)$ is a sinusoidal input. My understanding is that this system is described by the transfer function:
$$G(s) = \frac{x(s)}{u(s)} = \frac{1}{1+s/\omega_0}$$
I am interested in using this transfer function to calculate the amplitude and phase lag of $x(t)$. To do this, we evaluate $G$ at $s=i\omega$:
$$|G(i\omega)| = \frac{|x(t)|}{|u(t)|} = \frac{1}{\sqrt{1+\omega^2/\omega_0^2}}$$
$$\arg G(i\omega) = \text{phase lag} = -\tan^{-1}\frac{\omega}{\omega_0}$$
Are the above equations correct? Is my interpretation of them correct? Because when I numerically simulate the above differential equation in python, I get results that do not agree. Particularly, the blue dots and red dots should line up with the black lines in my plots. In other words, I have been unable to reproduce Fig. 3 in this tutorial

The ode is so simple, I don't see how I could have messed up the simulation. You can see my code in the links above. Thanks in advance for any comments/suggestions. (P.S. I hope I picked the right stackexchange site for this question, let me know if I got it wrong).

Update: I have tried decreasing the time step for the numerical integration routine, but this does not improve the fit. The code is still available at the same url: [Click here for code]

I think, the reason for the difference between both results could be the difference between the angular frequency $\omega=2\pi f$ and the "normal" frequency $f$.
That means: In the time domain you have defined the frequency $f$, but after processing you (perhaps?) have interpreted the results as $\omega$.


How to filter and downsample before FFT? I would like to perform STFT on musical signal with defined sampling rate (44100Hz).
I would like to get FFT result for smaller range of frequencies than are maximally available
(not 22050Hz but for 4000Hz) so the result will have better frequency resolution.
So I think I should downsample it, but how can you downsample to the frequency
that is not in integer relation to the input frequency (you cannot leave the kth sample)?
Do you do some kind of interpolation to get values of new samples?
How does it affect the signal?

I also know I should filter out the frequencies higher than 8000Hz otherwise
the signal would get aliased. Is there some kind of filter revelant for usage in STFT?

All I know about filters is that to get better (more sharp) filter
you need more coefficents and so it will take more time to calculate.
I know there are IIR filters and FIR filters.
I read IIR filters can be unstable but should I care about it if I use a ready implementation (and I assume it's done right)?

Is it better to filter the whole signal (whole audio file)
at once? If I get the infite response from the IIR filter and I filter
the whole signal at once will the energy be most smudged at the end part of the signal?

I also read filters can pose some kind of delay on the frequencies and it differs for different frequencies. How do you analyse and compensate this phenomenon?

Edit: I found that you can get FIR with linear delay and that now they are used more often than IIR. But there's still a delay, what does it mean? That the frequency events (like musical notes) will occur later in sample number time? What will be this delay?

Edit2: When I want to downsample to sampling rate of 8000Hz I have to filter so there won't be any frequencies over 4000Hz. Is this practically possible, because I looked at different FIR filters characteristics and they just seem to greatly damp the stopband frequencies not to eliminate them? If I will have just a bit of the high frequencies in the signal, will the signal look ok when downsampled (aliasing won't be noticeable)?

The process of changing the sampling rate by an arbitrary fraction $\frac{n}{k}$ is called resampling.  It is basically just interpolating and decimating at the same time.

It sounds like you are trying to go from 44100 Hz to 4000 Hz.  If so, that requires a sampling rate change of $\frac{40}{441}$, since $44100 * \frac{40}{441} = 4000$.  As a practical matter you don't want to do really large sample rate changes (like decimating by 441) all at once because it is very difficult to implement filters that can do that.  We thus try to break it down a bit into smaller steps.

If you factor the interpolation rate, 40, we see that it is 2 * 2 * 2 * 5.  We can likewise factor the decimation rate to 3 * 3 * 7 * 7.  One simple way to do the resampling, then, would be to do it in four steps like so:

$44100 Hz * \frac{2}{3} * \frac{2}{3} * \frac{2}{7} * \frac{5}{7}$

If you are trying to minimize the computational load you could rearrange the order of the resampling to reduce the sample rate as quickly as possible.  That would make the later resampling steps less computationally intensive because they don't have to work on as many samples.

$44100 Hz * \frac{2}{7} * \frac{2}{3} * \frac{2}{3} * \frac{5}{7}$

You could also combine the two $\frac{2}{3}$ resamples because their product, $\frac{4}{9}$, does not have any large numbers.  That leaves us with the following-

$44100 Hz * \frac{2}{7} * \frac{4}{9} * \frac{5}{7}$


How can I test an implementation of DFT for correctness? A friend of mine implemented a DFT algorithm for a project we're working on, and we're not quite sure how to test whether or not its correct. We're only in the eleventh grade, so we have very little experience in Fourier Analysis.

I remembered reading somewhere that the Fourier transform of the square wave function was the sinc function, so we thought that if the output of the program was the sinc function then the implementation was correct, but we quickly realized it wouldn't be that simple. Below is a picture of what the program outputs (blue is DFT scaled logarithmically to be in decibels):

!

Below is essentially the same as above only with lines instead of a solid fill:

!

Hmm, well that doesn't look anything like the sinc function. After looking at the code, I realized it was actually a plot of $\sqrt{r^2 + i^2}$ where $r$ and $i$ are the real and imaginary parts of the DFT respectively. So now I'm out of ideas as to how to test the validity of this program.

Does the output look correct? What other tests can we perform to verify the implementation?

A very quick and easy test that catches all relevant structural mistakes follows from a fundamental property of the Fourier transform (and also its inverse): It is a square root of the sequence order reversal operator. That means applying your implementation twice on a data vector will reverse the sequence order and applying it four times will give your original data back.

Note that this requires the so called unitary scaling convention, where the Fourier transform and its inverse share the same scaling factor. Otherwise you will get back a scaled version of the original signal.


Multi Scale Shape and Detail Enhancement from Multi Light Image Collections I have been working on a practical implementation in Matlab of the algorithm Raanan Fattal suggested in the following paper.

I couldn't understand the meaning of eq. 11 on page 5:

\begin{equation}\begin{array}{llllll}
U^{(i, j)}=g_{\sigma_{d}} * e^{\left(\left|D^{(i, j)}\right|-C^{(i, j)}\right)}.&&&&&(11)
\end{array}\end{equation}

g is a gaussian, 1D function (See section 4.1, page 3).
While in this equation no parameters are given.
$D(i,j)$ is a matrix and $C(i,j)$ is a number.

Could anyone explain what it means in case it is just a writing method I don't recognize?

Thanks

What the authors meant is to create a matrix of Weights, $ {U}^{\left( i, j \right)} $.

It is a matrix of the size of the image.
The given calculation is by the exponent of two terms (Each of them is a matrix, the calculation is element wise).
The final step is to blur it using a Gaussian Blue (2D Blur, like on an image).

In equation 10 you can see they take care of the normalization of the weights.

It would be nice if you come back and show us your results.


If MFCC vectors are speaker independent, then how can we reconstruct unique speech signals from them? Reconstructing a speech signal from a collection of MFCC vectors seems to work pretty well, but I've heard that one advantage of MFCCs is speaker-independence i.e. they are more-or-less the same across different speakers for a given phoneme. How then can a speech signal, with all its speaker-dependent idiosyncrasies (accent, etc.) be reconstructed from a supposedly speaker-independent MFCC vector? Are MFCCs not in fact speaker-independent then? If not, what does determine speaker-independence vs. speaker dependence?

Thanks.

First of all there is some serious "cheating" in the MFCC reconstruction experiment you linked to: not only the MFCCs are used, but also the voiced/unvoiced bit and the pitch.

MFCC are not speaker independent. In fact, they are used for speaker identification/verification tasks!

The speaker "idiosyncracies" are both in their prosody (preserved by this reconstruction experiment because the pitch is provided as a side-information to the reconstruction process) and in the articulation/timbre (preserved by the MFCC).

Two ingredients are needed to get MFCCs to work for speaker-independent recognition:

Vocal tract length normalization. A linear transformation (matrix multiplication of the MFCC vector) can map relatively well the MFCC sequences of two speakers speaking the same sentence. So even if MFCC are not speaker-independent one can optimize for a transformation matrix that "flattens out" the speaker-specific details.
Acoustic modelling. Using a large number of gaussians (or any classifier with large capacity) for a specific acoustic unit allows it to capture all the variations.


Determine the stability of a system without using the $\mathcal Z$-transform (described by a difference equation) For example, let's say a causal LTI System is described by the following equation:

$$y[n] - ay[n-1] = x[n] - bx[n-1],\quad n \in Z$$

Is there a way to determine (in this case) the stability of the system without using the $\mathcal Z$-transform?

(Solving the problem using the $\mathcal Z$-transform would be a classical easy problem - verify that ROC of $H(z)$ contains the unit circle, I know how to do that)

It boils down to the same thing, with or without $\mathcal{Z}$-transform. You determine the roots of the characteristic equation, which in your case is

$$\lambda=a\tag{1}$$

For stability, the magnitude of the roots of the characteristic equation must be smaller than $1$. So for your example you get the condition $|a|<1$ for stability.


How many multiplications and additions does it take to compute an FFT of a signal? How many multiplications/additions does it take to execute an FFT operation in terms of the length of the signal $N$ and the number of FFT bins $M$? I don't want an answer in terms of big O, but rather in terms of just $N$ and $M$.

If you are interested in more involved numbers of operations, you can check A modified split-radix FFT with fewer arithmetic operations, by Johnson and Frigo, IEEE Transactions on Signal Processing, 2007, authors of the "Fastest Fourier Transform in the West" (FFTW). For instance, they claim that the complexity of the classic radix-2 is about:
$$ 5N\log_2 N$$
while the split-radix yields:
$$4N \log_2 N − 6N + 8\,,$$
and they achieved the improvement:
$$\frac{34}{9} N \log_2 N − \frac{124 }{27 }N − 2 \log_2 N − \frac{2}{9}
(−1)^{\log_2 N} \log_2 N + \frac{16}{27}
 (−1)^{\log_2 N} + 8\,.$$
As said by others, you can see this number depends on implementations, and that the actual gain is subtle.


Does the FFT size matter in OFDM demodulation? The DFT size in OFDM demodulation is usually matched to the DFT size in the modulation process. This makes modulation/demodulation simple where each FFT bin (i.e., sample) represents one symbol. What is the impact of taking an arbitrary larger FFT size during the demodulation process (assuming properly sampled based on new DFT size)? Intuition seems to lead you to think that the underlying symbol is now spread across multiple bins. Can you use a larger FFT size and still demodulate the data?
EDIT: sample rate clarification

What is the impact of taking an arbitrary larger FFT size during the demodulation process?

It doesn't work.
This gets pretty intuitive if you take an arbitrary OFDM introduction and look at the spectrum of OFDM, consisting of overlapping sinc functions: you need to hit exactly that raster, or you don't see the individual channels, but an undisambiguable mishmash of carriers.
Using a larger FFT at the same sample rate can't work: don't forget that the FFT is just a DFT. A DFT maps vectors of a specific length N from time domain to frequency domain vectors of exactly the same size. So, a larger FFT length simply needs more time domain samples. But your OFDM symbol is as as many samples long as the transmit FFT was.
So, no, you can't use any other FFT size, sensibly.
You'd have to increase the sample rate by the same integer factor that you've increased the FFT size.
But that would only yield an FFT where all but the original bins are only noise, so you'd have won nothing.
(This is actually done in some OFDM systems to have more headroom for frequency synchronization, but it's not desirable if you don't have a frequency offset larger than your outer guard carriers.)


How to determine range of HSV values of the image? Is there any method to find range of hsv values of image (object) in rectangle inside big image? I mean, there is fixed size rectangle on my image(frame from camera) and i want to find range of hsv values of object which is inside this rectangle. Object is changeable. If i put red ball inside rectangle it should find range of red ball's hsv values, if i put my hand inside rectangle it should find range of my skin's hsv values. Is there any way to do this?

Edit:
I get part of image with ROI, and convert it to hsv format. But problem is behind red ball, I get wide variety of colours. So which range will I choose from this image?

Example: my frame:

Now I want to find best range of hsv values of image in rectangle box for tracking it. (Now it should be green and ranges should be h:53-90 s:74-147 v: 160-255)

If i can find this values, by using inRange method, i can convert it to binary image for tracking. Like this:

You can simply specify a ROI for that region and convert it into HSV.

For eg (below is pseudo-code in Python-OpenCV)

# define ROI of RGB image 'img'
roi = img[r1:r2, c1:c2]

# convert it into HSV
hsv = cv2.cvtColor(roi,cv2.COLOR_BGR2HSV)

Now it gives you the hsv values of the region.

But one or two difficulties there:

Your object may comprise some part of the ROI. There may be some background region also inside your ROI. ie if your object is red and background is blue, your ROI will contain two different values, one for RED and other for BLUE. Which one will your program select?
But colors won't be constant like RED or BLUE. You get wide variety of colours. Below is a plot of different colours for a ROI with green object.

But it has values from yellow to slight blue with green as major.

So you will have to select color with maximum probability, may be from checking its histogram.

Check this SOF for more details : https://stackoverflow.com/questions/13111504/selecting-best-range-of-values-from-histogram-curve


Filter order estimation Assume some unknown but small and finite number of poles and zeros in the complex Z plane, all with complex conjugates, producing some response.  Strictly from the absolute value of a set of equally spaced points around the unit circle, say greater than 2X the number of poles and zeros, of that response, is it possible to estimate or compute the number of poles and zeros which produced that sampled magnitude response?

Added:  Are more than 2X sample points required to determine the number of poles and zeros? (when given that the total is less than X).

Added:  If there is more than one solution, can a minimum solution (as in the minimum number of total poles and zeros) be found or estimated?

Theoretically it's possible to do this, although it will often not be practical.

Let's consider this in polynomial space. For a filter of order N you have 2*N+1 independent variables (N for the denominator and N+1 for the numerator). Let's look an arbitrary point $z_{k}$ in the z-plane and let's say the value of the transfer function at this point is H($z_{k}$). The relationship between the transfer function and all the filter coefficients can be written as equation that's linear in all filter coefficients as follows:
$$\sum_{n=0}^{2*N}b_{n}\cdot z_{k}^{-n} - H(z_{k})\cdot \sum_{n=1}^{2*N}a_{n}\cdot z_{k}^{-n}=H(z_{k})$$
So if you pick M different frequencies $z_{k}$ you'll end up with a set of M complex linear equations or 2*M real equations. Since your number of unknowns is odd (2*N+1) you probably always want to pick one frequency where z is real, i.e z = 1 or $\omega$ = 0.

If M is larger than N than the system of equations is linearly dependent. You can find the filter order by starting at N=1 and increase N until the equation system becomes linearly dependent. The largest N at which the system is linearly independent is the actual filter order. For this approach, it doesn't even matter what frequencies you pick. As long as they are different, any set of frequencies will work.

However, this a numerically very tricky problem. Polynomial representation for larger filter orders is numerically very fragile and the smallest amount of noise or uncertainty leads to very large numerical errors. For example, if you determine the values of the sampled transfer function through measurement, the required measurement accuracy will be prohibitive unless it's very benign low order filter.


FIR filter design with time and frequency domain constraints I deal with the problem of the synthesis of the filter with the requirements of both the frequency and
time characteristics.

Now there is a filter that meets the requirements.
Roughly speaking the requirements are as follows:

Frequency response

FIR Low Pass

Fpass ~ 0.0(3) (* pi rad/sample)
Fstop ~ 0.1(6) (* pi rad/sample)
bandwidth irregularity of ~ 1e-3 (maybe ~ 1e-2) dB
~ -40 dB suppression

Time requirements

Overshoot in step response ~ 4% (transition time ~ 1e-2 ms)

I would like to know about the possibilities of such a synthesis (I’m with filterDesigner in the matlab failed, tried free tools - Iowa Hills)
or at least about the transformation of the finished filter into a narrower
(Fpass = 0.025, Fstop = 0.125) (Matlab can do this, but only IIR is obtained).

The most straightforward way to impose constraints in the frequency domain as well as in the time domain for the design of (linear phase) FIR filters is to use linear programming.

Constraints on the step response or on the impulse response are naturally linear, and constraints on the frequency response can also be formulated as linear constraints in the case of linear phase FIR filters.

Take a look at this document for some examples. More details on the linear programming formulation of the FIR filter design problem can be found in this paper.


Where to get a proper dataset collection for bikes and cars I am doing a project on image classification and I don't know where to get the data-set collection for cars and bikes.Any suggestions to sites for this purpose is welcome.Also please suggest the general size of images to act as a database.

Thanks and regards in advance.

Bike and Cars database

http://www.emt.tugraz.at/~pinz/data/GRAZ_02/

http://homepages.inf.ed.ac.uk/rbf/CVonline/Imagedbase.htm

Hope this helps .. Download the bike.zip and cars.zip files and extract it.. If it is correct, click the tick mark to accept the answer.


How to reliably compute the group delay of a comb filter I applied the following FIR comb filter in real-time:

y[n]=x[n]-x[n-40]

Since this is an FIR, the group delay is D=(N-1)/2=20 samples. After applying the filter to a signal, I tried to use cross correlation between the filtered and unfiltered signal, to reproduce D computationally by determining the argmax of the cross correlation (I do have a need to the delay this way). The issue is that I get too peaks in the cross correlation, one at zero lag and another at 20 lag. But the peak at zero lag is the maxima which means the peak at 20 lag which is the correct lag is ignored. This method work really well with other filters like averaging filters.

Does anyone know while I get a the peak at zero which is overshadowing the real peak? Is this normal for comb filters? Is there another method to compute delays using the filtered and unfiltered signal other than cross correlation?

Since this is an FIR, the group delay is D=(N-1)/2=20 samples.

No, since this is a linear phase (i.e. symmetric or anti-symmetric) filter, the group delay is half the length! (being a FIR isn't sufficient.)

  The issue is that I get too peaks in the cross correlation, one at zero lag and another at 20 lag.

Write down the formula for auto-correlation at zero lag. Compare that to the formula of "energy of a signal". They are identical!

This really shouldn't surprise you!

  This method work really well with other filters like averaging filters.

This method works with anything that has a non-zero zero-lag coefficient.

  Does anyone know while I get a the peak at zero which is overshadowing the real peak?

Yes, because autocorrelation at zero is simply the energy. And since correlation is a linear, and your system passes through the original signal, plus a delayed version of it, you get the sum of the auto-correlation of the input signal and the cross-correlation of your delayed signal and the input signal.

The 20-lag peak is no "realer" than the 0-lag peak.

  Is this normal for comb filters?

This is normal for any linear time-invariant system.

  Is there another method to compute delays using the filtered and unfiltered signal other than cross correlation?

The group delay is really defined as derivative of the phase of your signal over frequency. If in doubt, estimate the spectrum of your system, and derive its phase. You'll notice that only a few specific systems (linear-phase, see above) have constant group delay.

Hence, I'm not sure your cross-correlation had much to do with group delay to begin with.


Initial value Problem with difference equation How to solve the initial value problem with

x(n + 1) + x(n) = n^2;
x(0) = 1

I do know how to solve for initial value theorem in differential equation, but not in difference equation.

Please help.Thanks.

Like a linear constant coefficient ODE the general solution can be written as the sum of a particular solution and the general solution of the homogeneous equation $x(n+1)+x(n)=0$.

Because the right hand side is a polynomial we look for a polynomial particular solution, and a bit of experimenting will show this to be a quadratic:$$x_p(n)=n^2/2-n/2$$

Also the general solution to the homogeneous equation is obviously: $x_h(n)=k(-1)^n$, and the solution to your problem is the sum $x_p(n)+x_h(n)$ with $k$ chosen to satisfy the initial condition.


eigen values and eigen vectors of signal What does the Eigen values and Eigen vector of a signal or function represent? What is its physical significance? I know about basis vectors of a signal which constitute the orthogonal planes where signal projections are represented. Are basis vectors and Eigen vectors same thing? Can we reconstruct signal using these Eigen vectors?

Consider a linear time-invariant system that maps a given signal to another signal space.  If the system produces a scaled version of the input signal $\phi$, say $\lambda \phi$, then we can view $\lambda$ and $\phi$ as eigenvalue and eigenvector respectively ($λ$ gives us the gain or attenuation of the eigen-signal).

Now suppose the impulse response of system is $h[n]$, when you input $x[n]$ is an eigen-signal, you have the output
$$y[n] = x[n]\sum\limits_{k=- \infty}^{\infty}  h[k]e^{-j \cdot \omega \cdot k}$$

so $$\lambda = \sum\limits_{k=- \infty}^{\infty}  h[k]e^{-j \cdot \omega \cdot k}$$

Note this is just the Discrete Time Fourier Transform of $h[n]$ since $H(e^{j \omega}) = \sum\limits_{k=- \infty}^{\infty}  h[k]e^{-j \cdot \omega \cdot k}$. Further, the Fourier Transform of $x[n]$ becomes meaningful as well.

Note that eigenvectors do not always form a basis. For example, $\begin{pmatrix} 0 &1 \\ 0 &0\end{pmatrix}$ has $0$ as its only eigenvalue, with eigenspace $\begin{pmatrix} x \\ 0 \end{pmatrix}$. There are not enough independent eigenvectors to form a basis.

For other discussions on the physical significance of eigenvalues or eigenvectors for a signal, please refer to this post from researchgate. And yes you can reconstruct the original signal using all the eigenvectors, or approximate the signal using some of them


Translation property of 2-D discrete Fourier transform I am studying the 2-D discrete Fourier transform related to image processing and I don't understand a step about the translation property.

In the book Digital Image Processing  (Rafael C. Gonzalez, Richard E. Woods    ) is written that the translation property is:

$$f(x,y)e^{j2\pi(\frac{u_0x}{M} + \frac{v_0y}{N})} \Leftrightarrow F(u-u_0, v-v_0)$$
and
$$f(x-x_0,y-y_0) \Leftrightarrow F(u,v)e^{-j2\pi(\frac{x_0u}{M} + \frac{y_0v}{N})}$$

So if I want to translate the Fourier Transform in the center of the center of the frequency rectangle, I had to consider $u_0=M/2$ and $v_0=N/2$.
So the first line becomes:

$$f(x,y)e^{j2\pi(\frac{M}{2}x\cdot\frac{1}{M} + \frac{N}{2}y\cdot\frac{1}{N})} \Leftrightarrow F(u-\frac{M}{2}, v-\frac{N}{2})$$

$$f(x,y)e^{j2\pi(\frac{x}{2} + \frac{y}{2})} \Leftrightarrow F(u-\frac{M}{2}, v-\frac{N}{2})$$

$$f(x,y)e^{j\pi(x+y)} \Leftrightarrow F(u-\frac{M}{2}, v-\frac{N}{2})$$

Then the book draws this conclusion:
$$f(x,y)(-1)^{x+y} \Leftrightarrow F(u-\frac{M}{2}, v-\frac{N}{2})$$

Why? What are the steps it did? I know that $j$ is a complex number then $j^2=-1$.

Very informally, they use the property $e^{j\pi} = -1$ (which is correct), the property $e^{ab}={(e^{a})}^b$ (valid for certain values of $b$), and derive the final expression for $b=x+y$ when $x+y$ is an integer.

The pattern $(-1)^{x+y}$ looks like a checkerboard.


Audio Broadband Noise Detection I have to quantify the broadband noise of an audio track.

I've seen lots of things around and the only viable option to estimate the broadband noise (without having any information about it) is something like that:

noise=0;
sig=0;
for i=in:fin
   nois(j)=(abs((min(B(:,i)).^2)));
   sign(j)=(abs((max(B(:,i)).^2)));
   j=j+1;
end
snr=var(sign)/var(noise)

Where B is the matrix of the spectrogram.

So, I'm doing FFT analysis of the audio, finding the minimum amplitude for each bin for an approximation of the noise floor and the maximum amplitude for an approximation of the signal. Then I'm doing the variances of both vectors and finally making a simple signal to noise calculation.

I know it's not a very sophisticated method but things previously tried have not worked. This gives me more solid results but it doesn't work for all kinds of genres.

Has someone something more to suggest me to improve this algorithm? Other approaches that work better?

The algorithm is basically based on the assumption that noise is "white", ie. in the expectation value, the noise energy in every DFT bin is equal. Under that assumption, and assuming a "long enough" observation, there has to be the same noise floor in every bin.

Sound on the other hand, is a highly correlated signal. So the assumption is that there's at least one DFT bin at a time where there's no significant sound signal.

As you might imagine, that's a bit of an oversimplification on both sides: There's a lot of music that actually uses Noise, or broadband signals, and audio noise typically is everything but white.

In essence: this boils down to the good ol' "What is my SNR?" question:

You, as the DSP designer, will have to model your signal, and your noise. After that, the question kind of answers itself. It's really more of a question of semantics than engineering: For example, take psychedelic rock/stoner/doom metal: wideband noise in the string instrument's track is very essential to that genre; is it noise, then? And if you get a recording of a classical guitar piece, wouldn't you consider the same noise in the signal part of the noise or music energy?


Is this a robust method for analysing correctness of a filtered signal? Is there a more robust method? I have a recorded signal which contains several frequency components with stochastic phase and amplitude changes and want to extract the true signal with the other frequency components and noise removed. I have performed several filtering techniques to attempt to achieve this. However I need a robust measure of how correct my filtered signal is and I do not have a method of measuring the signal on it's own, without the noise and other frequency components, in order to validate it that way. One method I have applied is minusing the filtered signal from the raw signal in time space and then taking the fft/psd (power spectral density) of the remaining signal and I see that the peak is gone, is this a robust measure of checking the signal I am tracking is correct? Is there any other methods for analysing this without having a way of measuring the true signal?

Some examples of the raw and differences signals for a couple of different filtering methods are shown below: (blue is the PSD of the original raw signal, and green is the PSD of [raw - filtered])

Given you can have a "reference waveform" representing the characteristics of your expected signals, to test the "correctness" of the signal extracted from noise use the Normalized Correlation Coefficient. I reference other posts below where I give the details of computing and using the Normalized Correlation Coefficient:

Noise detection

How can I find SNR, PEAQ, and ODG values by comparing two audios?

This will be an interesting approach to compare your different algorithms as the result of removing portions of your signal of interest or keeping noise contributions will both reduce the correlation coefficient as expected; so this will be an excellent approach to accurately indicate which algorithm you use is optimum for your waveform.


Matlab: Implement filter with no denominator A filter like this:
$z^{4}\cdot\frac{2 + z^{-1} - 3z^{-3} + ...}{1 - 0.9z^{-1}}$ seems hard for me to implement. The fraction part is easy:

filtered_signal = filter([2 1 -3 ...], [1 -0.94], signal);

But how is the $z^{4}$ part being implemented?

It's an odd notation but simple enough: $z^4% simeply means "shift the signal 4 samples to the left". That means its non-causal and can't therefore be implemented directly in a physical system.

In practice, the workaround that works most of the time is to add overall "causality delay" to the system until everything is causal again. Whether this is applicable here or not, depends on the specifics of your application.


What is the link between the homography computed over 2 images and the homography computed over the same images upside down? With OpenCV, I compute the homography between, say, these two images:

and

Don't worry about the strange white form on the right side, it is due to the smartphone holder I use. The homography, given by findHomography() function (using points detected with the Fast feature detector and the HammingLUT descriptor matcher), is:

A = [ 1.412817430564191,  0.0684947165270289,  -517.7751355800591;
     -0.002927297251810,  1.210310757993256,     39.56631316477566;
      0.000290600259844, -9.348301989015293e-05,  1]

Now, I use the same process to compute the homography between the same images that have been rotated by 180 degrees (upside down), using imagemagick (as a matter of fact, I would be equally interested to know the relation for rotation of 90 or 270 degrees...). Here they are:

and

With these images, the homography becomes:

B = [ 0.7148688519736168,  0.01978048500375845,  325.8330631554814;
     -0.1706219498833541,  0.8666521745094313,    64.72944905752504;
     -0.0002078857275647, -5.080048486810413e-05,  1]

Now, the question is how do you relate A and B? The two first diagonal values of A are close to the inverse of the same in B, but it's not very precise (.707805537 instead of 0.71486885). My ultimate purpose would to use the wanted relation to transform final matrix, avoiding to compute a costly image rotation.

I think when looking at the formula for the homography matrix:
$$
 H_{ba} = R - \frac{tn^T}{d}
$$
where $R$ is the rotation matrix by which $b$ is rotated in relation to $a$; $t$ is the translation vector from $a$ to $b$; $n$ and $d$ are the normal vector of the plane and the distance to the plane respectively (see Homography-3D plane to plane equation).

The relation between the two matrices lies in the normal vector of the plane. So you need to get the normal vector of the plane (out of the homography matrix), and apply the rotation to it, and then compute the homography matrix using the formula above. For the correct decomposing of the homography matrix, you can look at these code samples and this paper.


How to estimate filters using conjugate gradient? An image $I$ is computed by performing convolution and summation:

$$
\sum_{k=0}^{K-1} z_k * f_k = I
$$

Given only the feature maps $z_k$ and the resulting image $I$, how do I compute the filters $f_k$ ?

This is a linear system, therefore it must be possible to solve it with conjugate gradient. How to formulate it as $\mathbf{Ax} = \mathbf{b}$ ?

$z_k$ and $I$ have both the same quadratic size $N \times N$.

This is a typical deconvolution problem that you can solve either by transforming to the frequency domain where convolution is a simple multiplication:

\begin{equation}
F\left \{ {\bf{I}} \right \}= F\left \{ {\bf{z}} \right \} F\left \{ {\bf{f}} \right \},
\end{equation}

where ${\mathcal{F}}$ denotes the Fourier transform (or DFT) so $f$ will be:

\begin{equation}
{\bf{f}}={\mathcal{F}}^{-1}\left \{ \frac{{\mathcal{F}}\left \{ {\bf{I}} \right \}}{{\mathcal{F}}\left \{ {\bf{z}} \right \}} \right \}
\end{equation},
where ${\mathcal{F}}^{-1}$ is the inverse Fourier transform.

Otherwise, you can also represent in a matrix form, using a Circulant (kind of Toeplitz)  for the convolution operation. The frequency domain calculation though does not involve inverting (or pseudo inverting) a matrix.

If there's additive noise in your problem, please see Wiener deconvolution (both time and frequency domain implementation)


How do I know if my EEG signal need denoising? I recently started  working on sleep study.

For my research I download sleep EEG data from physionet. The EEG data has 100 Hz sampling rate and was recorded from 2 bipolar EEG site.

When I start the preprocessing stage, I encounter a simple problem:

How would I know if my signal has an artifact or noise?

It should be noted, based on Nyquist theorem and my signal's sampling rate, the maximum frequency of my signal is 50 Hz, so I did not filter unnecessary EEG frequency.

In general I only used a simple notch filter at 50 Hz, and used simple threshold method in order to remove the epochs that were grossly contaminated by muscle and/or eye movement artifacts.

Back to the main question, how should I know if I need to uses more complicated method for removing EMG or EOG artifact from my signal?

It primarily relies on,

what are you looking for in the EEG: define the frequency band that
relates with the phenomena you are interested to investigate and
filter out the rest (e.g. for sleep events bandpass ~= 0.8-30Hz).
power line of your area: notch filter at 50Hz or 60Hz.
Study design: if you are investigating events ensemble averages will remove random noise processes and highlight potential EEG events of interest.

Moreover, you can apply baseline removal to de-noise your signal from saturation trends and other systemical artefacts. Other experimental methods include ICA, which you will not be able to apply robustly since you only have two channels.