Autocorrelation for Stationary Signals I'm having trouble grasping the autocorrelation function for stationary signals, both strict stationary and WSS.

First for strict sense, we have $$\forall(\tau,t_1, \ldots, t_n) \in \mathbb{R} \land \forall n \in \mathbb{N} \ \ \ F_{X}(X_{t_1+\tau} ,\ldots, X_{t_n+\tau}) = F_{X}(X_{t_1},\ldots, X_{t_n}) \ \ $$ For unconditional cumulative distribution function $F_X$

This implies that all time instants are identically distributed $$\forall(t_1,t_2) \ \ F_X(X_{t_1})=F_X(X_{t_2})$$

My question then is, how is it possible for the autocorrelation to be anything other than a constant? The only thing I can think of is each time instant random variable $X_t$ is not independent, so the product distribution function would be different for different time instances and  $$E[X_{t_1}X_{t_2}]\ne E[X_{t_1}X_{t_3}]$$ if they were independent then the product distribution would be the same for all time instants and therefore the expected value as well giving a constant autocorrelation.

Additionally, for WSS, I can't conceptualize how autocorrelation could be independent of time, but still possibly dependent on shift. In this case however, each time instant is not necessarily identically distributed, only they have the same mean, but then it would seem the autocorrelation would depend on time.

Altogether, it seems that if the autocorrelation is time independent, then each time instant random variable is identically distributed (is this not the case?) and if they are identically distributed, the only way shift can affect the value is if

each time instant random variable has a different dependence on the other time instant random variables, but the dependence is consistent for all the random variables with respect to their "relative neighbors"

Is this understanding ok, or am I completely off the mark? Can someone provide an example where the autocorrelation of the stationary process is time independent, but does vary with shift?

You need the joint distribution of $X_{t_1}$ and $X_{t_2}$ in order to calculate $E[X_{t_1}X_{t_2}]$ and nobody has claimed that the joint distribution of $X_{t_1}$ and $X_{t_2}$ is the same as the joint distribution of $X_{t_1}$ and $X_{t_3}$. Stationarity requires that the joint distribution of $X_{t_1}$ and $X_{t_2}$ be the same as the joint distribution of $X_{t_1+\tau}$ and $X_{t_2+\tau}$ which is a quite different requirement than the one you need for claiming that the autocorrelation function of the process is a constant.  See this answer of mine for more than you might ever have wanted to know about autocorrelation functions.


Time shift -- two interpretations In their textbook(Example 1.1), Oppenheim and Willsky provide the following method for obtaining $g(t)=x(-t+1)$ from $x(t)$ : Replace $t$ with $-t$ in $x(t+1)$. In other words, advance the signal by 1 unit and flip the result.

But there could be another interpretation as well. Based on the fact that $-t+1 = -(t-1)$, we can have the desired signal as $h(t)=x(-t+1) = x(-(t-1))$,i.e. Delay the signal by 1 unit and flip the result.

But the two signals $g(t)$ and $h(t)$ are not identical. Why is the latter interpretation incorrect ?

Let

$$h_1(t) = x(t-1) ~~ \text{ for all }~ t, \tag{1}$$
that is, $h_1(t)$ is just
$x(t)$ delayed by one unit. Let $h(t)$ be the result of "flipping" $h_1(t)$, that is,
$$h(t) = h_1(-t) ~~ \text{ for all }~ t, \tag{2}.$$
But, $(1)$ says that $h_1(-t)$ is the same as $x(-t-1)$ for all $t$,
and so we conclude that
$h(t)$, the result of delaying $x(t)$ by one time unit and then flipping the
result gives us
$$h(t) = x(-t-1)~~ \text{ for all }~ t, \tag{3}$$
which of course is not the desired $x(-t+1)$.

So let's try O&W's prescription. Let $g_1(t) = x(t+1)$ for all $t$
be the result of
advancing $x(t)$ by one time unit. Now flip $g_1(t)$ to get
$$g(t) = g_1(-t) = x(-t+1)$$
which is what is desired, no?


characterization of DC component Consider the following two statements:

In time axis:

A signal without a DC component is a signal which doesn't have the zero frequency (the DC frequency)
A signal without a DC component is averaged to zero

My questions are:

Why is it incorrect? It sounds like the definition to me.
What is the explanation? It has something to do with the Fourier representation of the signal, I suspect.

Statement 1, as I read the garbled syntax, is correct, assuming "without a" and "doesn't have" imply a "zero" (an answer) and not a "null" (no answer).

The zero bin, also known as the DC bin, is the straight up sum of the signal over the sample frame in an unnormalized DFT calculation.  If the sum is zero, the average is zero.  So statement 2 is correct.

The explanation comes straight from the unnormalized definition of the DFT:

$$ X[k] = \sum_{n=0}^{N-1} { x[n] e^{ -i2 \pi { k \over N } n } } $$

Substitute in k = 0:

$$ X[0] = \sum_{n=0}^{N-1} { x[n] } $$

What is the source of these two statements?  Maybe there is a larger context to take into consideration.

Hope this helps,

Ced


Peak detection from saturated signals Say we have a sampled signal for which the inputs can be out of range. For instance, say a sampled sensor has a range of 0 to 1, and is presented with a time-varying stimulus where a peak occurs with max amplitude 1.2. In the data we would see a saturated or clipped peak.

What are some of the approaches to making estimates about the peak value?

How good of an estimate can be made by interpolating the value of the peak from the slope on either side of the clipping? This becomes difficult because if analysis relies on the dataset being at a specific sample rate, interpolation will result in upsampling, unless interpolation is done and the new signal is downsampled with an offset so that the peak becomes a data point. I'm sure this isn't an ideal approach.

What about calculations in the frequency domain?

Has anyone ever been presented with this problem?

Depends on the amount of clipping and what other information about the signal you have. It's pretty easy to detect clipping (wave form constant at +xmax or -xmax). Once a region has been found you can try to replace the clipped region with an interpolation from the non-clipped regions around this. The best interpolation scheme will depend on what else you know about the signal. You can try polynomial expansions, Taylor series, Frequency domain interpolation, etc.


How to change the start point of resampling process? I'm trying to resample a digital communication signal in SciPy or Matlab. Suppose the signal was sampled at 3 samples per second, and I want to resample it to 2 samples per second, the thing is I also want to be able to change the new sampling point so the new samples doesn't necessarily start at the same point as the orignal samples.

I'm using the Scipy resample function, but it doesn't have the flexibility to change the sample point, it only allows

The resampled signal starts at the same value as x ...

The Matlab resample function can't do that either.

Is there a way to shift the start point by a small amount, say $\alpha T$ where $T$ is the symbol duration and $0 \le \alpha < 1$ ? I guess what I need is the generic resampling algorithm underneath these resample functions so I can see if there's anything I can tweak.

since @Community resurrected this question, i'll pick up a little where hotpaw left off.  i'm not gonna run this through MATLAB or Octave (i can never remember which are columns and which are rows).

let $r$ be the resampling ratio.  $r>1$ means that your sample rate has increased and you're making the sequence longer because there are more samples.  so if you pad a delay of $\tau$ at the beginning of the input $x[n]$, that should result in a delay of $r \cdot \tau$ at the output $y[n]$ (whether $\tau$ is an integer number of samples or not).

now suppose you pad (using concatenation) the input with $P$ samples of zero before $x[0]$.  $P$ is a positive integer. that will result in a pad in the resampled output $y[n]$ of $rP$ samples where $rP$ is not necessarily an integer.  but

$$ \lfloor rP \rfloor \triangleq \operatorname{floor}(rP) $$

is an integer, and you can strip those samples (which will be mostly zero) off of the beginning of $y[n]$.  $\lfloor u \rfloor$ is the most positive integer that is no greater than $u$.  this is always the case:

$$ \lfloor u \rfloor \le u < \lfloor u \rfloor + 1 $$

at least for $u \ge 0$, we might say that $\lfloor u \rfloor$ is the "integer part" of $u$.  the "fractional part" of $u$ is whatever is left:

$$ 0 \le u - \lfloor u \rfloor < 1 $$

so, if what we want left is a fractional delay of $\alpha$ samples where $0 \le \alpha < 1$, then we want to choose an integer $P$ to pad the input such that $\alpha$ is the fractional part remaining in the output delay after stripping away the integer part, $\lfloor rP \rfloor$.  that is, given the sample rate ratio $r$, and fractional delay $\alpha$, we want to choose $P$ such that:

$$ rP = \lfloor rP \rfloor + \alpha $$

or

$$ \alpha = rP - \lfloor rP \rfloor \ . $$

after thinking about this a while, the only algorithm i can suggest is to start with $P=1$ and iteratively increment $P$ so that the above equation is close enough to be accurate.  that is, given an error constraint on $\alpha$, (call it $0 < \epsilon \ll 1$), increment $P$ until

$$ \left|\alpha - \left( rP - \lfloor rP \rfloor \right) \right| < \epsilon $$

pad the input with $P$ samples of zero, and strip the output of $\lfloor rP \rfloor$ samples, and you will be left with $\left( rP - \lfloor rP \rfloor \right)$ samples of delay.


How sampling aperiodic signal will result in periodic repetitions of the same I am reading "Digital Signal Processing" - Proakis and often read that sampling guarantees periodicity (not exact as read)

But I wonder how sampling aperiodic signal will result in periodic repetitions of the aperiodic signal?

sampling guarantees periodicity

Not exactly. Sampling in one domain guarantees periodicity in the other domain. So sampling in time creates a periodic spectrum and sampling in frequency creates a periodic time signal. Or the other one around: a periodic time signal has a discrete spectrum and a periodic spectrum has a discrete time signal.

That's why there are four different "types" of Fourier Transform. One for each permutations of continuous vs discrete in each domain. See for example http://fourier.eng.hmc.edu/e101/lectures/handout4/node3.html or https://www.dspguide.com/ch8/1.htm


May I know how to find Time frequency power spectrogram of a signal? I am analyzing EEG signal and I am unable to understand the term Time-Frequency Power Spectrogram. Please Can anyone help me?

Thank you in advance.

Short-time Fourier transform

Traditional Fourier analysis provides no information about evolution of frequency components over time. However, Short-time Fourier transform (STFT), a method which applies a short-time window to the signal and performs a series of Fourier transforms within this window as it slides across all the data, can overcome this limitation, providing a time–frequency representation of the signal.

In MATLAB, try

doc spectrogram

Continuous Wavelet Transform

Alternatively, continuous wavelet transforms (CWT) provide a useful approach in investigating non-stationary signals, this is often regarded as an “optimal” solution with regard to time and frequency resolution. Reportedly, CWT perform better at detecting abrupt changes in non-stationary signals compared to STFT.

In MATLAB, try

doc cwt

Generally, EEG signals have a low signal to noise ratio (SNR), therefore it is preferred to use ensemble averaging over trials for your time-frequency analysis (TFA), should the study design allow it.


Why calculate negative frequencies of DFT? This answer says that when calculating the DFT on a real valued signal, the (roughly) second half of the bins are the complex conjugates of the first half of the bins. Firstly I was curious if that is true?

If it is true, it seems like for frequency analysis purposes, you could forgo calculating the second half of the bins, which ought to be helpful for computation times and also for storage space if you want to store information in frequency domain since you could calculate the second half of the information on demand.

This also seems to be "more correct" to what people would expect when doing an FFT.

For instance, if using the equation:

$$X_k = \frac 1N \sum\limits_{n=0}^{N-1}x_ke^{-2i \pi kn/N}, \quad k\in[0,N), \quad k\in \mathbb Z$$

If you DFT a 4 sample cosine wave: $[1, 0, -1, 0]$, it gives the result: $[0, 0.5, 0, 0.5]$.

But, if you are only calculating the positive frequencies, it seems reasonable to modify the equation to be this:

$$X_k = \frac 2N \sum\limits_{n=0}^{N-1}x_ke^{-2i \pi kn/N}, \quad k\in[0,N/2), \quad k\in \mathbb Z$$

Which when applied to the 4 sample cosine wave gives us this result:
$[0,1]$

It seems more correct (to me anyways!) that it shows that $0\textrm{ Hz}$ (DC) has an amplitude of 0, and that $1\textrm{ Hz}$ has an amplitude of 1.  It's kind of confusing the other way, where the full amplitude of the $1\textrm{ Hz}$ wave is split between the positive and negative $1\textrm{ Hz}$ frequency.

Why is it then, that for real valued signals, we even bother calculating and reporting negative frequencies?

Mainly because its easier. The FFT is a specific algorithm to calculate the DFT. However, it only works if you calculate ALL frequencies (regardless if you want them or not). It takes in N complex values and it spits out N complex values. So the FFT can be used to evaluate the your first equation but not your second.

In your example, that's a trivial difference; but for larger FFT sizes this makes a substantial difference. It's a lot faster to calculate the values at negative frequencies with FFT and simply throw them away then calculating only the values at positive frequencies with a different algorithm.

This being said, there are ways for making use of the redundancy:

You can pack two real signals into a single complex, do a single complex FFT and split this into odd and even parts to get the DFTs of both input signals
There is a way to do a N point real DFT using N/2 point complex FFT. It just needs one extra stage of "unwinding" the result.


Determine impulse resonse of First Order Hold (FOH) Question, how can I determine the impulse response function of a first order hold? On Wikipedia it is simply stated as:

$$ h_{\mathrm{FOH}}(t)\,=  \frac{1}{h} \mathrm{tri} \left(\frac{t}{h} \right)  = \begin{cases} \frac{1}{h} \left( 1 - \frac{|t|}{h} \right) & \mbox{if } |t| < h  \\
0 & \mbox{otherwise}
\end{cases} \ $$

But given it's implementation
$$f(t) = f(kh) + \frac{t - kh}{h}(f((k+1)h) - f(kh)) \text{ with }kh \leq t < (k+1)h$$

how do I obtain $h_{FOH}(t)$?

--edit--

So if $f(kh)=δ(0)$ for $k=0$.

For $k=−1$, e.g., $−h≤t<0$ we get $f(t)=f(−h)+\frac{t+h}{h}(f(0)−f(−h))$ combined with $f(−h)=0$ this results in $f(t)=\frac{t+h}{h}\delta(0)$.

For $k=0$, e.g., $0≤t<h$ we get $f(t)=f(0)+\frac{t}{h}(f(h)−f(0))=δ(0)−\frac{t}{h}\delta(0)$.

So we get
$h_{\mathrm{FOH}}(t)\,=  \begin{cases} \frac{t + h}{h} \delta(0) & \mbox{if } -h \leq t < 0  \\ \delta(0) - \frac{t}{h} \delta(0) & \mbox{if } 0 \leq t < h  \\ 0 & \mbox{otherwise} \end{cases}$

what then follows, is that you have to "see" that this equals the formula with $\text{tri}$ function, as at the start of this message.

clear all;
close all;
clc;

h = 0.01;
t = -1:h:1;

for i = 1:length(t)
    if -h <= t(i) && t(i) < 0
        y(i) = (t(i) + h) / h;
    elseif 0 <= t(i) && t(i) < h
        y(i) = 1 - t(i)/h;
    else
        y(i) = 0;
    end
end

fohImpl = @(t,h) triangularPulse(t/h);

figure(1);
plot(t,y);
hold all;
plot(t,fohImpl(t,h));

Hint:

$f(t)$ is the expected output after passing a sampled function through the first order hold function with impulse response $h_{FOH}(t)$.

The sampled function is $f(kh)$, with $k$ being the sample index and $h$ being the sample interval, with zeros elsewhere as a continuous time process. The output is the continuous-time convolution of $f(kh)$ with $h_{FOH}(t)$.  In the first-order hold, $h$ is also the duration of the hold time, so spanning one sample interval.

If $f(kh) = \delta(0)$, with $k=0$, which is the unit impulse at time $t=0$, the result should be the impulse response, $h_{FOH}(t))$.


Making sense of the delay formula between two sensors in an antenna array? In the image above, there is a formula for the delay of a signal at a given sensor 'i' in an antenna array. I don't understand how they arrived at that formula.

Each antenna array element is separated by a distance of 'd' and we know the signal hits the array at an angle 'theta'. Dividing by 'c' gives us the time taken for light to travel the additional distance, but I don't see how the quantity 'dsin(theta)' represents the extra distance. Would appreciate any clarification.

but I don't see how the quantity 'dsin(theta)' represents the extra distance

Let analyse just sensors 1 and 2 also put $\theta$ 90 and 0 degrees and see what happens.

Let $\theta$ = 0 degrees, which means that transmitter of the signal is situated in front and both sensors should receive the signal at the same time (isn't it?).

$$
\Delta_i = { {d \cdot \sin {\theta} } \over {c} } = { {d \cdot 0 } \over {c} } = 0
$$

Let $\theta$ = 90 degrees, which means that transmitter of the signal is situated at the right hand and first sensors should receive the signal later. And delay should be time which need for signal travel the distance d with velocity c.

$$
\Delta_i = { {d \cdot \sin {\theta} } \over {c} } = { {d \cdot 1 } \over {c} } = {d \over c}
$$

But anyway in my opinion the indexes of sensors should be in the opposite direction.

I hope it helps you.


How should I modulate a $16$ bit input using FSK? I need to graph a $16$ bits binary sequence, assuming that it'll be a FM analog signal I need to apply FSK to it.

So, I think I should use the formula

$A \cos(2\pi f_1 t)$ for binary $1$'s

$A \cos(2\pi f_2 t)$ for binary $0$'s

But I'm really confused about how to set a method for getting the graph for my $16$ bit sequence ...

Since the comments on the question and the OP's responses
don't seem to be converging, here goes.

Pick a data rate $R$ bits/second (e.g. $10$ kilobits/second)
that is acceptable to your client or boss.
The FSK signal
will then convey one bit every $R^{-1} = T$ seconds (every
$T= 100$ microseconds if you chose the $R =10$ kbps data rate
suggested as an example), and the $16$ bits that you need to
transmit will be transmitted by an FSK signal $x(t)$ that will
have a total duration $16T$ seconds ($1600$ microseconds for the
example).  Specifically, let us assume that you need to define
what $x(t)$ is for all $t$ in the interval $[0, 16T]$.
Number the $16$ bits typed in by the user as $a_0, a_1, \ldots, a_{15}$
where each $a_i$ can have value $0$ or $1$. Include in your
program a snarky error message to be sent to the user if any of
the $a_i$ have value other than $0$ or $1$
Pick two positive integers $m_0$ and $m_1$. The FSK signal
described below will use frequencies $m_0/T$ Hz and $m_1/T$ Hz.
It is convenient, and possibly might even please your client or
boss, if you choose consecutive integers as $m_0$ and $m_1$ such
that the FSK signal fits in the frequency band allocated to the system.
Define the signal $x(t)$ as follows.  For $i = 0, 1, 2, \ldots, 15$,
$$x(t) = A \sin\left(2\pi\left[(1-a_i)m_0 + a_i m_1\left]\frac{t}{T}\right.\right.\right),~~~
iT \leq t < (i+1)T.$$
Several points need to be noticed with respect to this definition:

The quantity in square brackets has value $m_0$ if $a_i = 0$ and value $m_1$ if $a_i = 1$ and so during the time interval $[iT, (i+1)T)$, the signal
consists of either $m_0$ periods or $m_1$ periods of a sinusoid. At the
endpoints of the interval, the signal has value $0$. The signal thus has
the phase continuity that is desirable in an FSK signal as noted
by @JasonR in a comment.
If you are using some programming language to create an array of
$16$ gazillion points that you will then pass to a graphics package to display
as a graph of values of $x(t)$ as a function of $t$, then note that the
formulas used to create the first gazillion points will differ from the
formula used to create the next gazillion points etc.. The
formula needs to use $a_i$ to compute the values of the $i$-th
subset of gazillion points. If you want
to have a nice smooth graph of the sinusoids, it is necessary
to have $m_0$ and $m_1$ be onsiderably smaller than a gazillion
so that you get plenty of points
per period of each sinusoid.  Your graph will also have a long
horizontal axis that will likely extend past one page even if printed in
landscape mode on legal-size paper.


Finding the stability of the system Determine whether the system $h(n)=nu(n)$ is BIBO stable or not.

I have solved it and got it as unstable as $\displaystyle\sum_{n=-\infty}^\infty \lvert h(n)\rvert = \infty$. But the book ¹ is saying stable:

Which one is correct?

¹: S. Poornachandra, B. Sasikala: Digital Signal Processing. Third Ed. McGraw Hill, New Delhi

That's a (copy-paste) mistake in the book. The result in the given example is that the sum doesn't converge, yet the author concludes that the system is stable, which is obviously wrong:

$$\sum_{n=-\infty}^{\infty}|h[n]|=\sum_{n=0}^{\infty}n=\infty$$

The system is not BIBO-stable.


Calculate phase shift between two signals Consider the measurements of two modulated AC signals (voltage and current) which may have noise, harmonics and a DC bias. Each measured signal consists of a list for the time and a list for the corresponding values.

My goal is to find out the phase shift in degrees between these two signals in a robust way.

I would like to realize this in the frequency domain perhaps by calculating and evaluating the fundamental waves of the signals.

Any help or sample code (MATLAB) would be very much appreciated.

Do this procedure on both signals, and subtract the results to get what you want (it would be a nearly constant sequence whose value is phase shift in radians)

Perform FFT on the signal to obtain its spectrum.
Find the dominant harmonic of the signal : $f_c$ (where the peak of FFT occurs).
Filter the signal with a narrow-band linear phase filter (center = $f_c$, bandwidth as less as possible).
Again, perform FFT on the signal and extract its phase.

Note : It is assumed that the signal has a dominant frequency otherwise a phase shift is pointless


Unless I misunderstood completely, you are looking for the lag between the two waveforms red and green over smaller windowed intervals of that function as the window is moved over the sequence? In this case it can simply be done with a moving window using the following function in Matlab:

[corr, lag] = xcorr(a,b)

This will return the cross-correlation as corr, and the lag variable is the delay between a and b associated with each cross correlation.

Given you want to do this as a moving window, that can be achieved by setting a window size and index, and then index the window through your data:

Note this may be pseudo code below as I haven't tested the actual code written, but should give you the idea, as shown is a moving window of size K indexing 1 sample at a time, and will compute the relative delay between the waveforms within the window view.

x= green-mean(green);    % remove mean
x=x(:)';                 % ensure row vector
y = red-mean(red);
y=y(:)';
[~,n]=size(x);

for shift = 1:n-K    % ensure x and y same size
 xwin= x(shift:shift+K); % K is the window size
 ywin= y(shift:shift+K);
 [corr,lag]= xcorr(xwin,ywin);
 disp(shift)= lag(corr==max(corr));
end

If this is on track to what you are looking for, I would also consider normalizing the correlation by the standard deviation of the two waveforms within each window, thus giving you the correlation coefficient. This could allow for the option of optimally weighting the results for estimation of overall delay (or ignoring regions where the match was poor etc).


Is it possible to "equalise" a signal by deconvolving the impulse response of the room in which it is to be played? I am using a test sweep with a flat power spectrum and linear group delay (Optimized Aoshima's Time Stretched Pulse) to measure a room's frequency response. Having obtained the impulse response of the room, I thought I would be able to deconvolve it from my test signal, record the result in the same spot and under the same conditions, and obtain an approximately flat signal. However, the result of the deconvolution sounds like random noise, not at all like my original sweep.

I've tried naive deconvolution by inverse filtering, regularising my IR to avoid zeroes in certain frequencies, even ignoring all phase information in the IR and divinding just the magnitude of the results, all to no avail. In every case, the result of the deconvolution of sweep with the room IR sounds like random noise. Any ideas where I might've gone wrong?

Yes, the answer is here :
Compensating Loudspeaker frequency response in an audio signal ,
The filter here is your room transfer function.


Adding phase to a signal that has multiple frequencies An amplifier changes the phase at different amplitudes of asignal. I am using amplifier data in a MATLAB model.

I have constructed a signal that is a sum of sinusoids of different frequencies with random phase.

When this signal passes through the virtual amplifier model, how do I adjust the phase?

Should I... (Option 1)

Decompose the signal with a fast fourier transform, which gives me
the starting phase of each frequency and its absolute amplitude.
Then reconstruct the signal by making sinusoids of each frequency in
the signal (which will have the random phase it had to start), but
now add the additional phase change caused by the amplifier as well?

Or...  (Option 2)

Is this more complicated? If I have 3 frequencies in my signal, I
imagine 3 phasors (phase-vectors) added head to tail in a chain, with
the first one attached at an origin on a complex plane. They all
spinning around their tails at their own frequency. (so the first one
in the chain spins around the complex plane origin, the second one
spins around the head of the first, and the third spins around the
head of the second).
The absolute value and amplitude of my signal is
the vector from the origin to the end point of the final phasor in
the chain. Will the phase get added to this vector? Which is now what I said in the first part, the first option is where the phase is added to each phasor and not the overall phasor.

Please can someone assist

Caveat: many behavioral models of amplifiers exist, none of which is perfect, and I don't know what is best for your particular amplifier.  But the most common model is basically your option 2.  I have never heard of anyone using option 1.

Let $x(t)$ be the input to the amplifier in complex baseband form.  A very common amplifier model is to assume that the gain is a memoryless function $g(|x(t)|)$ of the input amplitude, and the phase shift is also a memoryless function $\phi(|x(t)|)$ of the input amplitude.  The complex baseband output is therefore

$y(t) = g(|x(t)|)\exp(j\phi(|x(t)|))x(t) $

This model assumes the signal is relatively narrowband.

There is no need to decompose the signal into its Fourier coefficients to perform this transformation, and indeed this only complicates the matter.


Help regarding property of unit impulse function We know that using properties of unit impulse function it can be shown that

$$\int_{t_1}^{t_2} x(t) \delta^{(n)}(t-t_0) dt=(-1)^nx^{(n)}(t_0),\quad t_1<t_0<t_2$$ (source: Continuous and Discrete Signals and Systems- Soliman)

But what would be the case for the problems like these
$$\int_{-1}^{2}(3t^3+9)\delta''(t+1) dt$$

I've tried in this way:

$$\int_{-1}^{2}x(t)\delta''(t+1) dt=\\ \big[x(-1)\delta'(t+1)-x'(-1)\delta(t+1)\big]_{-1}^{2}-\int_{-1}^{2}x'(t)\delta'(t+1)dt$$

can't figure out the first part inside the third bracket.

Integrals of the form

$$\int_{a}^bf(t)\delta^{(n)}(t-a)dt\tag{1}$$

with $\delta^{(n)}(t)$ being the $n^{th}$ generalized derivative of the Dirac delta impulse, are undefined. What is well-defined are integrals of the form

$$\int_{a^+}^bf(t)\delta^{(n)}(t-a)dt\tag{2}$$

or

$$\int_{a^-}^bf(t)\delta^{(n)}(t-a)dt\tag{3}$$

In $(2)$, the Dirac impulse (or its derivative) is outside the integration interval and, consequently, the integral evaluates to zero, whereas in $(3)$ $\delta^{(n)}(t-a)$ is inside the integration interval and the value of the integral is given by the first formula of the question.


How can I specify the distance between values in such row I have a question, about verifying the value which are near to each other in such row, Is there any algorithm which can determine that ?

Suppose I have a vector $z_i =$ {$z_1 , z_2, z_3, . . . ., z_n$}, I want first to determine the maximum and minimum values of this vector, then check if the distance between each other value is smaller compare then maximum and to minimum values. In other words, it means  dividing the vector into two categories and then check if the value is nearer to first category or second.

Is there any algorithm can do that?

thank you

Just a little bit of python:

import numpy as np

def splitVector(v):
    thres = (np.max(v)-np.min(v))/2
    return v[v<=thres], v[v>thres]

This will yield two vectors with just the lower or higher elements. If you want to preserve overall size and position of arguments, just a little modification will do it:

import numpy as np

def splitVector(v):
    thres = (np.max(v)-np.min(v))/2
    return np.where(v<=thres,v,0), np.where(v>thres,v,0)

This will fill in zeros for all elements not fulfilling the condition.


Pixel color labeling tool Anyone know good tool for manual pixel color labeling like SegNet do?

LableMe does same work as SegNet but manually. You can learn about this tool from this video tutorial

To start using this tool you have to create an account on LabelMe website.


What is meant by a system's "impulse response" and "frequency response?" Can anyone state the difference between frequency response and impulse response in simple English?

The impulse response and frequency response are two attributes that are useful for characterizing linear time-invariant (LTI) systems. They provide two different ways of calculating what an LTI system's output will be for a given input signal. A continuous-time LTI system is usually illustrated like this:

In general, the system $H$ maps its input signal $x(t)$ to a corresponding output signal $y(t)$. There are many types of LTI systems that can have apply very different transformations to the signals that pass through them. But, they all share two key characteristics:

The system is linear, so it obeys the principle of superposition. Stated simply, if you linearly combine two signals and input them to the system, the output is the same linear combination of what the outputs would have been had the signals been passed through individually. That is, if $x_1(t)$ maps to an output of $y_1(t)$ and $x_2(t)$ maps to an output of $y_2(t)$, then for all values of $a_1$ and $a_2$,

$$
H\{a_1 x_1(t) + a_2 x_2(t)\} = a_1 y_1(t) + a_2 y_2(t)
$$

The system is time-invariant, so its characteristics do not change with time. If you add a delay to the input signal, then you simply add the same delay to the output. For an input signal $x(t)$ that maps to an output signal $y(t)$, then for all values of $\tau$,

$$
H\{x(t - \tau)\} = y(t - \tau)
$$

Discrete-time LTI systems have the same properties; the notation is different because of the discrete-versus-continuous difference, but they are a lot alike. These characteristics allow the operation of the system to be straightforwardly characterized using its impulse and frequency responses. They provide two perspectives on the system that can be used in different contexts.

Impulse Response:

The impulse that is referred to in the term impulse response is generally a short-duration time-domain signal. For continuous-time systems, this is the Dirac delta function $\delta(t)$, while for discrete-time systems, the Kronecker delta function $\delta[n]$ is typically used. A system's impulse response (often annotated as $h(t)$ for continuous-time systems or $h[n]$ for discrete-time systems) is defined as the output signal that results when an impulse is applied to the system input.

Why is this useful? It allows us to predict what the system's output will look like in the time domain. Remember the linearity and time-invariance properties mentioned above? If we can decompose the system's input signal into a sum of a bunch of components, then the output is equal to the sum of the system outputs for each of those components. What if we could decompose our input signal into a sum of scaled and time-shifted impulses? Then, the output would be equal to the sum of copies of the impulse response, scaled and time-shifted in the same way.

For discrete-time systems, this is possible, because you can write any signal $x[n]$ as a sum of scaled and time-shifted Kronecker delta functions:

$$
x[n] = \sum_{k=0}^{\infty} x[k] \delta[n - k]
$$

Each term in the sum is an impulse scaled by the value of $x[n]$ at that time instant. What would we get if we passed $x[n]$ through an LTI system to yield $y[n]$? Simple: each scaled and time-delayed impulse that we put in yields a scaled and time-delayed copy of the impulse response at the output. That is:

$$
y[n] = \sum_{k=0}^{\infty} x[k] h[n-k]
$$

where $h[n]$ is the system's impulse response. The above equation is the convolution theorem for discrete-time LTI systems. That is, for any signal $x[n]$ that is input to an LTI system, the system's output $y[n]$ is equal to the discrete convolution of the input signal and the system's impulse response.

For continuous-time systems, the above straightforward decomposition isn't possible in a strict mathematical sense (the Dirac delta has zero width and infinite height), but at an engineering level, it's an approximate, intuitive way of looking at the problem. A similar convolution theorem holds for these systems:

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

where, again, $h(t)$ is the system's impulse response. There are a number of ways of deriving this relationship (I think you could make a similar argument as above by claiming that Dirac delta functions at all time shifts make up an orthogonal basis for the $L^2$ Hilbert space, noting that you can use the delta function's sifting property to project any function in $L^2$ onto that basis, therefore allowing you to express system outputs in terms of the outputs associated with the basis (i.e. time-shifted impulse responses), but I'm not a licensed mathematician, so I'll leave that aside). One method that relies only upon the aforementioned LTI system properties is shown here.

In summary: For both discrete- and continuous-time systems, the impulse response is useful because it allows us to calculate the output of these systems for any input signal; the output is simply the input signal convolved with the impulse response function.

Frequency response:

An LTI system's frequency response provides a similar function: it allows you to calculate the effect that a system will have on an input signal, except those effects are illustrated in the frequency domain. Recall the definition of the Fourier transform:

$$
X(f) = \int_{-\infty}^{\infty} x(t) e^{-j 2 \pi ft} dt
$$

More importantly for the sake of this illustration, look at its inverse:

$$
x(t) = \int_{-\infty}^{\infty} X(f) e^{j 2 \pi ft} df
$$

In essence, this relation tells us that any time-domain signal $x(t)$ can be broken up into a linear combination of many complex exponential functions at varying frequencies (there is an analogous relationship for discrete-time signals called the discrete-time Fourier transform; I only treat the continuous-time case below for simplicity). For a time-domain signal $x(t)$, the Fourier transform yields a corresponding function $X(f)$ that specifies, for each frequency $f$, the scaling factor to apply to the complex exponential at frequency $f$ in the aforementioned linear combination. These scaling factors are, in general, complex numbers. One way of looking at complex numbers is in amplitude/phase format, that is:

$$
X(f) = A(f) e^{j \phi(f)}
$$

Looking at it this way, then, $x(t)$ can be written as a linear combination of many complex exponential functions, each scaled in amplitude by the function $A(f)$ and shifted in phase by the function $\phi(f)$. This lines up well with the LTI system properties that we discussed previously; if we can decompose our input signal $x(t)$ into a linear combination of a bunch of complex exponential functions, then we can write the output of the system as the same linear combination of the system response to those complex exponential functions.

Here's where it gets better: exponential functions are the eigenfunctions of linear time-invariant systems. The idea is, similar to eigenvectors in linear algebra, if you put an exponential function into an LTI system, you get the same exponential function out, scaled by a (generally complex) value. This has the effect of changing the amplitude and phase of the exponential function that you put in.

This is immensely useful when combined with the Fourier-transform-based decomposition discussed above. As we said before, we can write any signal $x(t)$ as a linear combination of many complex exponential functions at varying frequencies. If we pass $x(t)$ into an LTI system, then (because those exponentials are eigenfunctions of the system), the output contains complex exponentials at the same frequencies, only scaled in amplitude and shifted in phase. These effects on the exponentials' amplitudes and phases, as a function of frequency, is the system's frequency response. That is, for an input signal with Fourier transform $X(f)$ passed into system $H$ to yield an output with a Fourier transform $Y(f)$,

$$
Y(f) = H(f) X(f) = A(f) e^{j \phi(f)} X(f)
$$

In summary: So, if we know a system's frequency response $H(f)$ and the Fourier transform of the signal that we put into it $X(f)$, then it is straightforward to calculate the Fourier transform of the system's output; it is merely the product of the frequency response and the input signal's transform. For each complex exponential frequency that is present in the spectrum $X(f)$, the system has the effect of scaling that exponential in amplitude by $A(f)$ and shifting the exponential in phase by $\phi(f)$ radians.

Bringing them together:

An LTI system's impulse response and frequency response are intimately related. The frequency response is simply the Fourier transform of the system's impulse response (to see why this relation holds, see the answers to this other question). So, for a continuous-time system:

$$
H(f) = \int_{-\infty}^{\infty} h(t) e^{-j 2 \pi ft} dt
$$

So, given either a system's impulse response or its frequency response, you can calculate the other. Either one is sufficient to fully characterize the behavior of the system; the impulse response is useful when operating in the time domain and the frequency response is useful when analyzing behavior in the frequency domain.


Wavelet filter coefficients from the Morlet wavelet How can I obtain the wavelet filter coefficients ( i.e.       lowpass decomposition filter,        highpass decomposition filter,        lowpass reconstruction filter,        highpass reconstruction filter) for the Morlet wavelet?

For convenience: the wavelet curve follows the equation $$\pi^{1/4}  \exp(jwx) \exp(-0.5x^2)\,.$$

The critically-sampled  lowpass/highpass and subsampled filter bank implementation of the Discrete Wavelet Transform (DWT) is only valid for a certain class of wavelets. And the Morlet wavelet family does not belong to this class.

Explained differently, discrete wavelets seldom have closed-form expressions, because of the constrains laid on the dyadic scales and filters. Since the Morlet wavelet is complex, you may want to invest in one of its discrete avatars, the dual-tree wavelet transforms, that are somewhat redundant however. I happen some personal references on the topic:

Noise Covariance Properties in Dual-Tree Wavelet Decompositions
Image Analysis Using a Dual-Tree M-Band Wavelet Transform
A nonlinear Stein based estimator for multichannel image denoising


Do all real-time filters introduce phase lag, or can there be such a thing as "lag-less" or "zero lag" real time filters? Do all real-time filters introduce phase lag? I once heard of "lagless filters", but I don't think they can be used real-time (ie: on live data streaming in for a control system, for example, operating on incoming data), correct?

Obviously even though this question could be simply answered in a binary "yes/no" fashion, expounding upon the answer will be very helpful.

Update:

I just added some links in the comments below this question.

Also, here's some brief comments from a couple controls experts I know:

Kenny Jensen:
"Forward-backward filters, aka smoothers, have no lag. You can't run them online though."

Ed Lim: "Hence they are not "real-time" filters. So I believe the answer is yes regarding real-time filters."

If a filter has any memory and that memory affects the current output, then that past data is lagged into the present.  If that filter is time invariant, then a memory of the present can potentially affect future output. e.g. is lagged into the future.  So, yes, unless the filter is memoryless or not time invariant.

Note that a unit impulse at the origin is memoryless, but impossible to implement in any finite size system due to the speed of light being finite.


Notation confusion -- What is the correct operator for computation of the log-likelihood expression for complex valued data? This question is an extension of another question of mine asked earlier here Help in understanding if the maximum likelhood estimation is working properly

In that question the inputs were real valued. Now, in this question the inputs are real valued and hence the confusion.

I am having difficulties in the notation when applying to the complex domain.
What notations and operators should I use when all the data are complex valued? Please help

Without seeing your overall signal model we cannot really help. I suppose $y$ is your estimate? And hence, $L=L(y)$. Where is $h$ in your Likelihood equation? However, some general observations:

As Batman pointed out: likelihood corresponds to probability, which is necessarily a real value. So, if your equations yield complex values, you did something wrong. (Note that, depending on numerical accuracies etc. it can happen, that there is a very small imaginary part ($\approx 10^{-16}$) in your numeric computations. This part can be ignored.)
Commonly (but not always!) when in the real-domain you have a transpose in your equation, it changes to conjugate transpose (i.e. Hermitian operation) in the complex domain.
Noise variances commonly remain real-valued. However, noise realizations will be complex-valued.
You should understand the expression in your exponential as the squared distance between the measurement and the estimate. In real domain, the squared distance of $x,y$ is equal to $(x-y)^2$. However, in the complex domain, it becomes $\|x-y\|^2=(x-y)(x-y)^*$.


How to Make a Gaussian Filter? The 2D Gaussian function is defined as:
$$G(x,y) = \frac{1}{2\pi \alpha}e^{-\frac{x^{2}+y^{2}}{2\alpha}}$$
And this is a Gaussian kernel:

If I take the top left corner as the origin and set $\alpha=1$, then at $x$=4 and $y$=0, $G(4,0)$ = $5.3\times10^{-3}$. But here at $G(4,0) = 3.66\times10^{-3}$. Then how they design this Gaussian kernel? I am new on this topic. So, this question may seem too much beginner level.

Cool facts about the Gaussian surface:
It is a rotation:
$$ G(x,y) = \frac{1}{2\pi \alpha}e^{-\frac{x^{2}+y^{2}}{2\alpha}} = \frac{1}{2\pi \alpha}e^{-\frac{r^{2}}{2\alpha}} = G(r) $$
where $ r = \sqrt{x^2 + y^2} $
It is separable:
$$ G(x,y) = \frac{1}{2\pi \alpha}e^{-\frac{x^{2}+y^{2}}{2\alpha}} = \frac{1}{2\pi \alpha} e^{-\frac{x^2}{2\alpha} } e^{-\frac{y^2}{2\alpha} } = G_x(x)G_y(y) $$
With proper normalizations
$$ G_x(t) = G_y(t) = G_t(t) = \frac{1}{\sqrt{2\pi \alpha}}e^{-\frac{t^{2}}{2\alpha}} $$
This means every slice of a Guassian surface is a Guassian function.

I used to do a lot of smoothing on scatter dot diagrams to make them nice surfaces.
The first is the same as DC.  You want the sum of your entries to equal the denominator you are using.  Yours is fudged a little bit from the nearest solution found here.
Here is the closest match for a denominatory of 271.  You can see where they tweaked to get 273.

[[  1.   4.   6.   4.   1.]
 [  4.  16.  26.  16.   4.]
 [  6.  26.  43.  26.   6.]
 [  4.  16.  26.  16.   4.]
 [  1.   4.   6.   4.   1.]]

[[  0.78997163   3.54040722   5.83714469   3.54040722   0.78997163]
 [  3.54040722  15.86700434  26.16026756  15.86700434   3.54040722]
 [  5.83714469  26.16026756  43.13098958  26.16026756   5.83714469]
 [  3.54040722  15.86700434  26.16026756  15.86700434   3.54040722]
 [  0.78997163   3.54040722   5.83714469   3.54040722   0.78997163]]

The bottom one is the actual surface rescaled by the denominator.
I am wondering like Royi in the comment if this is supposed to be an exercise in integer approximation.  So, I got curious and without further comment on the issues that are raised by this program, here it is:

import numpy as np

#==========================================================
def main():

        N = 5

        alpha = 1.0

        theSurface = np.zeros( ( N, N ) )

        L = ( N - 1 ) / 2

        C = 1.0 / ( 2.0 * np.pi * alpha )

        for x in range( N ):
          dx  = x - L
          dx2 = dx * dx
          for y in range( N ):
            dy  = y - L
            dy2 = dy * dy
            r2  = dx2 + dy2

            s = C * np.exp( -r2 / ( 2.0 * alpha ) )

            theSurface[x,y] = s

        print( theSurface )

        theIntegerMatches = []

        for theDenominator in range( 100, 10000 ):
          theApprox = np.round( theDenominator * theSurface )

          theSum = np.sum( theApprox )
          if abs( theSum - theDenominator ) < 0.5:
             theIntegerMatches.append( ( theDenominator, theApprox ) )

        for theMatch in theIntegerMatches:
          theDenominator = theMatch[0]
          theApprox      = theMatch[1]

          theModel = theApprox / theDenominator

          theDiff  = theSurface - theModel
          theDiff2 = np.multiply( theDiff, theDiff )

          theDiffSum  = np.sum( theDiff )
          theDiff2Sum = np.sum( theDiff2 )

          print( theDenominator, theApprox[0][0], theDiffSum, theDiff2Sum, 1.0 / theDiffSum )

        print

        theMatch = theIntegerMatches[-1]
        theDenominator = theMatch[0]
        theApprox      = theMatch[1]

        theActual = theDenominator * theSurface

        print( theApprox )
        print( theActual )

#==========================================================
main()

Any formula based on r^2 is going to be a rotation around an axis.
That being said.  In integer implementations it is customary to choose a power of two denominator so that the division can be performed with a simple bit shift.
Keep this in mind when dealing with any filter like this:  Your filter is going to work on a surface that can be approximated by the following Taylor series representation.  At any point, there is a best fit set of coefficients.  In this case, your filter zeros out several of these.
$$
\begin{aligned}
f &= C_1
\begin{bmatrix}
1 & 1 & 1 & 1 & 1 \\
1 & 1 & 1 & 1 & 1 \\
1 & 1 & 1 & 1 & 1 \\
1 & 1 & 1 & 1 & 1 \\
1 & 1 & 1 & 1 & 1 \\
\end{bmatrix}
+ C_{x}
\begin{bmatrix}
-2 & 1 & 0 & 1 & 2 \\
-2 & 1 & 0 & 1 & 2 \\
-2 & 1 & 0 & 1 & 2 \\
-2 & 1 & 0 & 1 & 2 \\
-2 & 1 & 0 & 1 & 2 \\
\end{bmatrix}
+ C_{y}
\begin{bmatrix}
-2 & -2 & -2 & -2 & -2 \\
-1 & -1 & -1 & -1 & -1 \\
 0 &  0 &  0 &  0 &  0 \\
 1 &  1 &  1 &  1 &  1 \\
 2 &  2 &  2 &  2 &  2 \\
\end{bmatrix}
\\
&+ C_{xx}
\begin{bmatrix}
4 & 1 & 0 & 1 & 4 \\
4 & 1 & 0 & 1 & 4 \\
4 & 1 & 0 & 1 & 4 \\
4 & 1 & 0 & 1 & 4 \\
4 & 1 & 0 & 1 & 4 \\
\end{bmatrix}
+ C_{yy}
\begin{bmatrix}
4 & 4 & 4 & 4 & 4 \\
1 & 1 & 1 & 1 & 1 \\
0 & 0 & 0 & 0 & 0 \\
1 & 1 & 1 & 1 & 1 \\
4 & 4 & 4 & 4 & 4 \\
\end{bmatrix}
+ C_{xy}
\begin{bmatrix}
 4 &  2 &  0 & -2 & -4 \\
 2 &  1 &  0 & -1 & -2 \\
 0 &  0 &  0 &  0 &  0 \\
-2 & -1 &  0 &  1 &  2 \\
-4 & -2 &  0 &  2 &  4 \\
\end{bmatrix}
\end{aligned}
$$
These are not the only basis matrices, but they are probably the easiest to understand.


Memorylessness of simple delay system As usual, $y(t)$ is the output signal of a system, and $x(t)$ is the input signal.  I'm wondering whether or not a certain system has memory.

It's easy for me to see that the system

$y(t) = \int_{t-T}^{t} x(t) dt$ where $T > 0$

has memory, because an output at time $t_0$ depends on all values of input from $[t_0 - T, t_0]$.

But what about the following system, which just delays input by some constant $C$:

$y(t) = x(t - C)$ where $ C > 0 $

Upon first glance, this doen't feel like it has memory because there is no integral.  But the output of a delay system does depend upon inputs from the past.  Does that mean a simple constant delay system has memory?

A system is memoryless if its output ($y(t)$) for each value of the independent variable ($t$ in this case) at a given time is dependent only on the input at that same time ($x(t)$).

Every system that consists of a delay (like the one in your example) or an accumulator, for example, are systems with memory.

This can be seen just by replacing $t$ by some number. Take, for instance, $t=0$. In that case, the output for that instant ($y(0)$) depends on the input for another value of the independent variable ($x(C),C>0$); that means that the system has memory.


Measuring how "non-noisy" a sound signal is Given a sound clip (whose source is an FM radio signal) I'd like to measure how different from background noise (the "in-between channels" noise) the sound clip is, so that I can compare the "non-noisiness" (and thus the quality) of two radio signals. Is there a good way to do this?

If the noise is (close to) stationary, you could use one of the many noise-floor estimation methods to get an estimate of the noise power and, consequently, of the SNR. Search for "noise floor estimation" and "noise floor tracking". You'll find a lot of publications and patents.


CWT at low scales: PyWavelets vs Scipy Low scales are arguably the most challenging to implement due to limitations in discretized representations. Detailed comparison here; the principal difference is in how the two handle wavelets at each scale:

pywt: diff of int_psi resampled at each scale, self-normalized (caveat2); int_psi == cumsum(psi)
scipy: wavelet recomputed and explicitly normalized at each scale

(For pywt it's not exactly diff(int_psi), but key behavior is equivalent) Below plots compare real and complex morlet CWT's for a sharp transition signal; all coefficients were L1-normalized.

I didn't get much from inspecting shapes of respective wavelets (blue = recomputed, orange = resampled, both rescaled -1 to 1):

So what's the deal? Am I missing something in the discretization differences, or is pywt's integrated wavelet approach superior for lower scales (if so, why)?

Verdict: I conclude scipy's higher leftmost peak is due to both pywt's wavelets' lower amplitude at lower scales and scipy's wavelets' stronger correlation with lower frequencies at lower scales. Overall, they both correlate approximately evenly with irrelevant frequencies (scipy being a little worse), but pywt "lucks out" with its bad correlations being offset with smaller wavelet amplitudes at lower scales, thus the smaller peak.
I'm not exactly sure why scipy's direct recomputation is worse than pywt's resampling; I'd say resolution problem with very limited number of samples early on, and resampling integrated working out better (it's also slightly longer), but haven't investigated explicitly.
For additional assessment on low scales pywt vs scipy, see bottom section. The rest of the answer serves to justify the verdict, but doesn't add much further; the most informative is Fig. 2. Some code is provided, which can be reused to test other transforms.

A more naive kind of norm points to the answer; zero-mean + scale coeffs between $[-1,1]$ per scale:

For the top plot, a $[0,1]$-rescale is equivalent, so it shows the normalized correlation strength with the wavelet irrespective of energy (actually not quite; caveat1). This suggests the difference to be rooted in wavelet amplitudes; let's compare. (-- coeffs code3)

Do they correlate differently? -- For all plots, blue = scipy, orange = pywt.
Scale-standardized: divide each scale vector by own maximum (highest absolute correlation among all timeshifts). Plot absolute value of each scale vector, so timeshifts along horizontal. A greater value at a given timeshift indicates greater correlation of the plotted scale with relative to all other timeshifts at that scale.

Above shows scipy's early scales having comparable correlations at each scale for both halves of the signal, whereas for pywt the right half is dominant.
Does this imply scipy correlates more strongly with lower frequencies at lower scales than pywt (i.e. for same wavelet amplitude)? Yes; we'd get identical scale-standardization if we forced wavelet amplitudes the same. See zoomed amplitude-equalized plot in Fig 2; note how at indices 2 and 4 pywt exceeds scipy, and observe the above gif.

Timeshift-standardized: divide each timeshift vector by own maximum (highest absolute correlation among all scales). Plot absolute value of each scale vector, so
timeshifts along horizontal. A greater value at a given timeshift indicates greater correlation of the plotted scale relative to all other scales at the timeshift.

Above shows, for left half, about same correlation strength for all 20 scales relative to strongest-correlating scale. Right-half fluctuates much more for both, but scipy keeps itself multiples above pywt. This implies most of the energy discrepancy at these scales is due to the right-half.
But why due to right-half? And why is left-half relatively fixed? Figure 1; toward higher scales (lower frequencies, left-half), scipy and pywt agree in both wavelet amplitude and shape.

Wavelet amplitudes comparison
Instead of looking at max amplitude, I define a measure of "mean amplitude": mean of absolute value of tail-trimmed wavelet, where "tail" = any absval 1e7 times less than peak amplitude (instead of strictly zero which is rarer). This is to unbias the mean for wavelets with long tails: (-- code2)
left  = [i for i, v in enumerate(x)       if abs(v) > x.max() / 1e7][0]
right = [i for i, v in enumerate(x[::-1]) if abs(v) > x.max() / 1e7][0]
meanamp = np.mean(np.abs(x[left:-(right + 1)]))

At scales[0], the ratio is 2.58 - when squaring for the energies, that's a x6.64 larger coefficient for scipy due to wavelet amplitudes difference alone. However, the ratio of coefficients' energies at that scale is 8.33, so there's another factor - and this factor is, surely, scipy's stronger correlation with low frequency part of the signal, as shown earlier.
To see how much exactly the amplitude accounts for, we can normalize by it; scale coefficients by ratio of (mean) amplitudes:

We see for most lower scales scipy's energies greater even with amplitude accounted for, leaving only correlation strength as the source of difference. Also interesting to see pywt correlate much more strongly at the first true mode (~scale 50).

Caveat1: Fig. 0 is misleading
Notice scipy jump above pywt at the right edge, yet, if we look at logscale mean amplitudes, pywt's is both greater, and correlates stronger as per the amplitude-normed plot. What's the deal? Recall how scale-standardized plot was defined; this is same, but now for coefficients - scale vectors are rescaled relative to themselves, rather than between pywt and scipy.
A greater value in Fig 0 results if most points aren't far below the maximum (as all points are divided by max) across timeshifts, or that the wavelet at the scale correlates more consistently for each timeshift. So it's more of a consistency measure; here's the last scale, shown for every timeshift:

So... pywt's better at lower scales? -- Uncertain, but pywt wins in two important test cases:
Higher max frequency

Noisy, combined frequency signal


How can I detect the tree rings? I'm looking for the way to detect the tree rings automatically in the wood dist image.
Using the matlab, I have removed the background, and binarized the image.

How can I detect the white and black boundaries?
I want to fill that white area.

Is there a way to detect the white and black boundaries?
I want to get rid of the black noise that exists in the white area.
is there any other way to make clearer binary image?
It has scratches which has caused by sandpaper machine, so morphological operation doesn't work well.. any other way..?

I would appreciate your help.

It looks like you could use a low pass filter or another noise reduction filter. Then maybe an edge detection filter.


Why is White Noise so important in System Identification or Adaptive Filters I'm looking to implement a feedback cancellation filter using Wiener Filtering, where an adaptive Wiener filter is used to cancel the feedback occurring in the path between a loudspeaker and a mic (assume PA system). The idea is essentially from this paper:

  Spriet, Ann, et al. "Adaptive feedback cancellation in hearing aids
  with linear prediction of the desired signal." IEEE Transactions on
  signal processing 53.10 (2005): 3749-3763.

According to the paper:

Signal delivered at the mic is $x[k]$.
Signal driving the loudspeaker is $u[k]$.
The transfer function of the feedback path is $F[q]$, so the signal captured at the microphone is $y[k] = x[k] + u[k]F[q]$.
The feedforward transfer function from microphone to the loudspeaker is $G[q]$.

Here is the block diagram that helps put it all together:

In the description, $k$ is the discrete time index, and $q^{-1}$ is the unit delay operator (I know these to conventionally be $n$ and $z^{-1}$).

The idea is to introduce an adaptive filter $\hat{F}[q]$ that estimates the feedback path and cancels it from $y[k]$. There's a couple of things being talked about, but essentially since $x[k]$ and $u[k]$ are correlated, they talk of adding a probe signal $r[k]$ to the loudspeaker input $u[k]$, which helps identify $\hat{F}[q]$. $r[k]$ is apparently usually a noise, which I assume helps decorrelate $x[k]$ and $u[k]$.

  My first problem is: Wouldn't the noise also be output by the
  loudspeaker? Or is it added to a copy of the signal, not affecting what is fed to the loudspeaker?

Alternately, the paper also says that many audio signals can be closely approximated as a low order AR process:

$x[n] = h[n]*w[n]$, where $w[n]$ is white noise. This condition is not satisfied for voiced speech or music, in which case a pulse train is suggested.

  So my bigger question is, why is this noise (white noise) so important
  in AR modeling, or Adaptive Filtering? It seems to defeat the purpose
  to add noise to the signal.

Any

I believe the point of feeding white noise into the system is for the filter to adapt its coefficients before actually generating the signal $x[k]$. This would mean there are two "operating modes" for the system: coefficient adapting mode (in which white noise, a broadband signal, is used to adapt the filter to the feedback path), and performing mode (where the signal $x[k]$ is input into the system with the filter already adapted).

So the answer to your first question would be that the noise is output by the loudspeaker and follows the whole signal path. The system adapts it's coefficients to be similar to $F(q)$, and once that process has ended the white noise can be turned off. Afterwards, when generating $x[k]$ the system will already reduce feedback using the $\hat{F}_0(q)$ transfer function that was adapted previously, and the coefficients will only change when the forward path $G(q)$ is modified.

The purpose of using white noise (or a pulse train) is to have a broadband signal that has frequency components throughout the whole audible spectrum so that the adapted filter works for any kind of input signal $x[k]$. If the system was used without first feeding the broadband signal into it, every time $x[k]$ presents new frequency components, the coefficients may need to adapt during the performance/speech, which may have audible effects depending on the convergence speed of the algorithm.


Interpolation of audio for new frames I am trying to upsample a video by existing frame interpolation techniques. In the process, I realize that I also need to interpolate the audio signal for the new frames so that the audio signal is in sync with the frames of the new video.

The new upsampled video should have the same frame rate as the original input video, and will be a slowed down version of the input video.

Is there any way I can proceed to do the part of the audio interpolation?
If the frame rate is F and the audio sampling rate is Fs, then Fs/F audio samples are associated with each frame of the original video.

Any insight on how to proceed regarding this idea would be helpful.

When upsampling the number video frames to allow playing a video in slow motion, rather than for a frame rate increase for smoothness, you will likely need to modify the audio using a time-pitch stretching/shifting algorithm to stretch the audio out (increase the number of samples) for a longer play duration (to match the increase in slow motion video playback time), and to do so without distorting vocal or musical pitches.

If you just interpolate or resample audio without time-pitch correction, you will mess up all the pitch frequencies in the audio.

See this Wikipedia article for a description of several possible time-pitch stretching/scaling algorithms. There are multiple commercial and open source applications and libraries that can accomplish time-pitch modification/correction with varying degrees of quality and artifacts.


Dominant eigenvectors of an unknown matrix Do you have any idea about how we can find the principal eigenvectors of an unknown matrix ${H}$?

The elements of $H$ are unknown in general. If you are familiar with channel estimation procedure in wireless communications, $H$ is channel matrix, and every element $(i,j)$ represents the channel gain between antenna element $i$ of transmitter and antenna element $j$ of the receiver. The following are the information that we have:

we do not know the channel matrix,
we know second order statistics, covariance matrix of the channel,
the channel matrix size tends to infinity (large number of antennas in Tx and Rx),
the channel is sparse, that is, it has only two or three dominant eigenvectors regardless of the channel matrix size.

Note that the power method is not a good solution for me, since its complexity is huge, especially because it requires many iterations to find the dominant eigenvectors.

Some information about the interplay between channel and its covariance:

We consider Rayleigh correlated channel coefficients $h_k\sim CN (0,\mathbf{R}_k)$, mutually independent across the users $k$, and $CN$ is circularly symmetric complex Gaussian. We let $\mathbf{R}_k = \mathbf{U}_k \mathbf{\Lambda}_k \mathbf{U}_{k}^{H}$, where $\mathbf{U}_{k}$ is a tall matrix of the eigenvectors corresponding to the non-zero eigenvalues of $\mathbf{R}_k$, given as the diagonal elements of the diagonal matrix $\mathbf{\Lambda}_{k}$.

The Karhunen-Loeve representation of the user channel vector is given as $\mathbf{h}_{k} = \mathbf{U}_{k}\mathbf{\Lambda}_{k}^{\frac{1}{2}}\mathbf{w}_{k}$, where the entries of $\mathbf{w}_{k}$ are i.i.d. with zero mean Gaussian distributed with variance 1.

In the absence of other answers, and sorry if you have already considered this, but the power method can be very fast with a sparse matrix. You may like to check out how Google does it (eigenvectors are used for page rank, which is pretty fast).

http://en.m.wikipedia.org/wiki/PageRank

  Because of the large eigengap of the modified adjacency matrix
  above,[22] the values of the PageRank eigenvector can be approximated
  to within a high degree of accuracy within only a few iterations.

Footnote 22 points to: http: //www-cs-students.stanford.edu/~taherh/papers/secondeigenvalue.pdf

Footnote 23 discusses how to make the calculation even faster
http: //projecteuclid.org/euclid.im/1150474883

I have not read the papers, but I hope that you find one of them to be useful. (I had to break the links as new users can not post many links.)


How to determine the filter type of a discrete frequency spectrum I have the discrete impulse response of a filter and i want to determine which kind of filter it is by using the DFT.

The impulse response is:

h0 = [0,1,2,1,0]

which I zero-pad for a better result

N0 = 10;
h1 = [zeros(1,N0),h0, zeros(1, N0)];

Then i created the N-Point DFT of it using FFT

H1 = fft(h1);

The picture shows the plot of abs(H1)

Now i know that k1 = 0 represents the steady component and at first I thought that this is the spectrum of a band-stop filter but now that I read a little bit more about the DFT I'm not quite sure because the DFT is always mirrored for real signals.

So my question is: Do I have to use only one half of the DFT spectrum to determine the filter characteristics? (Which would mean that this is a low-pass filter?)

The $k_1$ indexes you get from the FFT correspond to frequencies 0 to $2\pi$ (actually, 1 sample short of $2\pi$). Note that frequencies in the range $\pi$ to $2\pi$ also correspond to the negative half of the specturm ($-\pi$ to 0).

You can get a nice visualisation with frequency 0 in the middle by applying fftshift to the result (Octave or Matlab assumed).

You can also check freqz, which plots a nice frequency and phase repsonse (try freqz([1 2 1], [1])). Or try plot(abs(freqz([1 2 1],[1]))) to get a linear plot rather than log.

What you are getting is a low-pass filter (gain 4 at low frequency, 0 at high frequency).

In VERY general terms: You can know beforehand it is a LPF because the impulse response is a single bump. A high-pass filter will have some sort of antisymetric shape (like [1 2 -2 -1]), and a band-pass or band-reject will have some oscillations in the impulse response.


Implementation of modulation/demodulation models I am trying to implement mod demod for data rate of few kbps . I would like to know is it better/ easier to implement in DSP/FPGA. I am aware of both C and VHDL programming. I think the code flow and build and test time would be faster in DSP? please suggest? Any particular resource for implementation (like a sample code) to get started?

I would suggest you to use the GnuRadio API usualy used for Software Defined Radio (SDR) projects.
As Marcus Müller said in the comments, it is easier to develop this kind of project and test it with a software implementation.
I am also pretty sure that you will find a lot of existing examples to help you with your project.
Guided Tutorials
API documentation (not beginner-friendly)
Source code


How do I create a frequency vs time plot? I'm a chemical engineer, not an EE, so this is a bit difficult.

I'm trying to figure out how to take amplitude vs time data and transform it into frequency vs time.  My first instinct is to slice my data into chunks, perform the FFT on each chunk, and then plot that.  Unfortunately, as the time duration of each slice approaches zero, there is no longer enough information to get accurate frequency information (low frequencies require more than a very small time slice).  So... how do I do this?  I'm sure this is some sort of famous problem that someone's already solved.

Here is the kind of transform I'm looking for, illustrated with a sound wave (piano note G).  As you can see, this graph as three axes, with the third one being represented by color.

Thanks!

The time vs frequency resolution is a well-known problem, and there are indeed approaches to overcome it. For audio signals, some of the commonly used techniques include: parametric methods ; adaptive resolution (analyze with various time/frequency configurations and patch the results together - Wen X. and M. Sandler, "Composite spectrogram using multiple fourier transforms") ; wavelets/decompositions on overcomplete bases ; and use of phase information to extract the precise location of frequency peaks (IFgram).

However, it appears that the graph you have shown does not use some of these techniques ; so I suspect this is not what you might be looking for. There appears to be some "smearing" on the horizontal axis (for example at t=1.2s) and this is a sure sign that the analysis has been done with a high overlap between chunks.

Indeed, the chunk duration and the number of analysis frames per second do not have to be linked to each other if you allow frames to overlap. So if you want to use 40ms long analysis frame, your grid does not have to be:

frame 1: t=0..t=40ms ; frame 2: t=40ms..t=80ms

It could very well be:

frame 1: t=0..t=40ms ; frame 2: t=10ms..t=50ms

This overlapping can give the illusion of a higher temporal resolution without reducing too much the FFT window size. Note that this can only help in accurately locating an event on the time axis - it will not help resolving two events close in time... Just like increasing the FFT size might help with identifying the location of a frequency peak, but not with the resolution of two adjacent frequency peaks.


Downsampling - how many extra samples to include I want to downsample a finite section of an equidistant time-series
$$x_i, \qquad i=0\dots n-1$$
 by doing an FFT and dropping the upper part of the spectrum. Since the FFT looks at $x_i$ as a cyclic series, it will show edge effects (ringing) near $i=0$ and $i=n-1$ caused by for example a jump between $x_{n-1}$ and $x_0$ in the cyclic view. Hence I need to include some extra samples left of $x_0$ and right of $x_{n-1}$, i.e. take the indices
$$i=-m\dots n-1+m$$
into the FFT, where the corresponding samples $x_{i<0}$ and $x_{i>n-1}$ are either accessible values (measurements) or respresent a reasonable extrapolation (e.g. zero-padding, $x=0$). The extra samples are dropped again once the downsampling is completed.

Now the question is, how many samples should I add? It is clear that the 'sinc' function, which represents the antialiasing cutoff (and hence, the edge ringing) is infinitely extended. So, I can, of course, not expect to obtain an exact solution from adding a finite set of samples. On the other hand, some quick examples show that adding just one sample on either side is not enough in most cases, since this only goes one half swing of the 'sinc' function towards the boundary values.

What heuristics could lead me to choose a reasonably appropriate number of extra samples (probably by specifying an accuracy parameter)?

The image above shows the result of a downsampled time-series that originally started at 1, but has been extended by a bunch of zeroes to the left. Since the ringing is not only determined by the sample rate but also by the height of the original step (in this case =1), it is not so easy to determine when the ringing (sinc) has dropped 'enough'. Imagine the step scaled by 1000, then the ringing is also scaled by 1000.

I want to downsample a finite section of an equidistant time-series
  by doing an FFT and dropping the upper part of the spectrum.

Why? This one of the least efficient way to down sampling.  See Why is it a bad idea to filter by zeroing out FFT bins?

In order to down-sample you need to low pass filter first. There are many different ways of building a low pass filter and almost all of them are better than zero-ing FFT bins.

  What heuristics could lead me to choose a reasonably appropriate number of extra samples (probably by specifying an accuracy parameter)?

You do this by first specifying the relevant requirements for your application: signal to noise ratio, phase distortion, pass band ripple, stop band ripple, latency, CPU/memory constraints. Some of those may be very important, some of those may be irrelevant, but it's different for each application.  Then you design a filter that makes sure you meet your requirements.


Justification for Squared $ {L}_{2} $ Data and Smoothness Term as an Error Bound Often in variational methods (and not only) we have an energy that is of the form:
$$E(u) = \frac{1}{2}\|f-u\|^2_2 + \frac{\alpha}{2}\|\psi(u)\|^2_2,$$

where the first term is referred to as the data term, and the second as the smoothness term. I understand that the squared $L_2$ norm is especially appealing since it results in a simple expression when we are looking for the minimum. I am looking for the formal motivation for designing such an energy. More specifically, I would like to see it as resulting from some bound upon the actual error. Let $f$ be a degraded/noisy image of the original image $g$, and $u$ be the image that we are looking for. Is it possible for the above to be interpreted as a bound of $\|g-f\|^2_2$? Let us consider an even simpler problem (where $h$ is a sufficiently smoothed version of $f$):
$$E(u) = \frac{1}{2}\|f-u\|^2_2 + \frac{\alpha}{2}\|u-h\|^2_2$$
If I start from $\|g-f\|_2$ I can get the following bound:
$$\|u-g\| = \|u-f+f-g\| \leq \|u-f\| + \|f-g\|$$
$$\|u-g\| = \|u-h+h-g\| \leq \|u-h\| + \|h-g\|$$
Combining the two inequalities with some weight $\lambda \in [0,1]$ I get:
$$\|u-g\|\leq (1-\lambda)\|u-f\| + \lambda\|u-h\| + (1-\lambda)\|f-g\| + \lambda\|h-g\|$$
Since $f,g,h$ are fixed, the minimization is only over the first two terms:
$$\min_u (1-\lambda)\|u-f\| + \lambda\|u-h\|$$
And this can be rewritten through $\alpha$ by setting $\alpha = \frac{\lambda}{1-\lambda}$. My issue is that these are still not squared norms. If I try to square both sides, this results in a different minimization energy. The other way to get the squares involves writing out something of the form:
$$\|u-g\|^2 = \|(u-f) + (f-g)\|^2 = \|u-f\|^2 + \|f-g\|^2 - 2(u-f)\cdot(g-f) \\ \leq \|u-f\|^2 + \|f-g\|^2$$
The last inequality holds only when $(u-f)\cdot(g-f) \geq 0$ though. Said otherwise, the angle between $(u-f)$ and $(g-f)$ needs to be smaller than 90 degrees which is not something that can be guaranteed I believe.
A looser bound is:
$$\|u-g\|^2 \leq \|u-f\|^2 + \|f-g\|^2 + 2\|u-f\|\|f-g\|,$$
but it also involves non-squared terms.

Here's something else I noticed about the non-squared energy:

$$\alpha\|u-f\| + \beta\|u-h\|, \alpha, \beta > 0$$

will always have as solutions $u=f$ or $u=h$ for $\alpha<\beta$ respectively $\alpha>\beta$. When $\alpha=\beta$ any solution of the form $(1-\lambda)f + \lambda h$ is admissible. It is thus clear that the above is not very interesting as a minimization energy.

Cross-posted to: https://math.stackexchange.com/q/3692290/463794
I expect to get a different perspective there.

One of the motivations to use the $ {L}_{2} $ norm comes from the Maximum a Posteriori Estimation (MAP) framework.

If you model $ \psi \left( u \right) \sim \mathcal{N} \left( 0, \alpha \right) $ then if you derive the MAP Estimator in case the added noise is Gaussian you'd get the exact model you posted above.

An example of the derivation of MAP model to the above can be seen in my answer to Estimating the Signal by Deconvolution with a Prior on the Filter Coefficients and the Signal Samples.


Programming a parametric equalizer - preventing clipping I've been developing a music player for Android for many years, and I've recently started working on the implementation of a parametric equalizer.

I'm having trouble preventing a gain at one band from adversely affecting the gain at a neighbouring band, resulting in a peak gain higher than desired, which ultimately may result in clipping.

I've been following the Orfanidis paper on 'Digital Parametric Equalizer Design With Prescribed
Nyquist-Frequency Gain' as well as an example Matlab implementation of some of the concepts outlined in the Orfanidis paper.

The equalizer effectively acts as a graphic eq - I've fixed the center frequencies and bandwidths to specific values:

(Values in Hz, Bandwidth in brackets)

63 (15)
125 (30)
250 (30)
400 (45)
630 (75)
1000 (115)
1600 (180)
2500 (290)
4000 (460)
6300 (730)
10000 (115)
16000 (1825)

The plan is to allow a maximum boost/cut of 12dB, with bandwidth gain set to 3dB below peak. Ideally, I could first apply a 12dB cut across the entire spectrum, so that a peak boost of 12dB is guaranteed not to cause any clipping.

The problem is, if a user applies a boost of 12dB at, say, center frequencies 63Hz and 125Hz, the peak gain at both frequencies ends up being more than 12dB, which may result in clipping. Here's what the FFT looks like for that case:

Notice the peak gain is actually greater than 12dB.

The effect is compounded if a boost is applied across all the frequencies mentioned above (this is the worst case scenario):

Notice the peak gain is more like 15dB.

Maybe I should expect this, and the answer is simply that bands shouldn't be spaced too closely if applying a lot of gain. But, my problem is I'm trying to provide this functionality to a casual music listener, but I don't want them to configure the equalizer in such a way that clipping occurs, as they'll likely claim that it's faulty.

So my questions are - does this look like a correct implementation of a parametric eq? Does it make sense that neighbouring bands can have a greater peak gain than the specified 'boost gain'? How can I prevent a user from causing clipping in this situation?

Any help is appreciated.

There's no universally acclaimed correct design of equalizers. Just like there are multiple car designs per brand, there are also multiple techniques for achieving desired effects.

It's true that your channel filters will interact with each other and change the effective resulting gain at each frequency, but mostly pronounced at the adjacent channel peak though.

The amount of interaction depends on your selected DC gains (which are unity for your application) and channel bandwidths and bandwidth gains.

There are recommended practices to prevent distortion levels. Headrooms are included so that peak gains will not create clipping. Or you can allow some...?

Robert Bristow Johnson is a member here and he's one of the designers of those (cookbook) parametric equalizers that you want to use in your application.


How can I automatically identify and remove a fundamental frequency and all of its harmonics? I have some microscopy data that is contaminated by a heartbeat artifact that I'd like to remove. The data consists of a large time series of images captured at ~60Hz.
Here's a small example clip in GIF form:

I have taken the average pixel intensity over time, and computed the periodogram using Welch's method:

As you can see there is a sharp peak at ~1.8Hz which is likely to correspond to the heart rate (~108 beats/min). There are also a bunch of harmonic peaks at integer multiples of 1.8Hz. The exact heart rate is likely to vary from dataset to dataset, but I can specify a biologically plausible range as shown by the shaded area on the periodogram.
What I'd like to be able to do is:

Automatically detect the fundamental frequency corresponding to the heartbeat, and all of its harmonics
Filter the data so as to remove the fundamental and all harmonics.

At the moment I can solve point 1 very crudely by finding the largest peak in the periodogram, then multiplying it by $1, 2, ..., N$ where $N$ is the estimated number of harmonic peaks, but I'm sure that there must be a better method than this hack.
Regarding point 2, I came across this question which mentions using a comb filter to remove a fundamental and all of its harmonics. Is this the best method to use? One important consideration is that I will have to apply the filter to each pixel timeseries in a large array, so a computationally efficient method would be highly desirable.
Example data

I've uploaded a 1D timeseries of example mean
pixel values in .csv format here.
I've also added a sample of 1000 frames (spatially decimated 2:1) in a compressed .mat file here

Your method isn't that bad for a first attempt.

However, the following method tends to work better:

Search for local maxima
Check for maxima that are close (2 or 3 bin spacing) and merge them
Create a few hypotheses about the fundamental frequency. You currently assume the highest peak is the fundamental frequency, which is one hypothesis. You should also check for the possibility that the highest peak is the first harmonic, i.e. that there's a smaller peak at half the frequency of the main peak. You may also have other cases to consider, using your knowledge of the problem at hand (Interference? Irregular heartbeat?).
Assuming each of those hypothesis, find the ground frequency by fitting parabola to each harmonic peak. Each peak will produce a slightly different estimate due to noise, but these errors are uncorrelated and average out. One of the hypotheses will lead to much better fits, pick the ground frequency predicted by this hypothesis.
Using the ground frequency you found in step 4 as a given, refit the parabola around each peak to estimate the height of the peak. Note that the peaks probably fall between two bins.
You now have the location and strength of the fundamental and its harmonics, but not the phase. It's probably easiest to find the phase of the fundamental, subtract that, find the phase of the first harmonic, etc.

The core reason that this works better is step 4. Any wrong hypothesis for a ground frequency will fail horribly as you're trying to fit peaks around the location of predicted harmonics. Say you've got a peak at 2 Hz. This could be the ground frequency or a first harmonic. When you test the "first harmonic hypothesis", i.e. whether the ground frequency is actually 1 Hz, you fit parabola's to the data around 1,2,3,4,5 ... Hz. If that hypothesis is wrong, you get garbage around 1,3,5 Hz. If it's right, you might find peaks near 1.1 Hz, 2.2, 3.3, 4.4 and 5.6 - which suggests that the actual ground frequency is 1.12 Hz.


what is Shannon entropy? I was reading a paper about finding optimal morlet wavelet function for CWT, and they were minimizing Shannon entropy.
what is Shannon entropy and it uses?
(specially on continuous wavelet transform)

Shannon entropy is a way of measuring the degree of unexpectedness or unpredictability of a random variable. For example rolling a die has higher entropy than flipping a coin because the die has more possible outcomes making it harder to predict. Same goes for a biased coin versus a fair coin.


Image processing research areas using GPU I'm an undergraduate student and I am currently trying to find suitable topic for my bachelor thesis. The broader topic I've selected is "Image processing on GPU", so basically I'll be developing image processing algorithms specifically for GPU.

The problem is, I need to do a novel approach, not reinvent a wheel.  There's so much going on in image processing, lately on GPUs as well, that I'm having trouble finding a suitable area that hasn't been researched in this way and which isn't above my capabilities as an undergraduate.

The question is, can you point me in a direction what are current research areas in image processing that could benefit from GPU implementation? I'd love to do something for medical or astrophysical applications, but I'm open to other areas as well.

Thanks!

I've done a bit of GPGPU (for optical flow) during my PhD, so here are a few hints to help you:

choose a computation intensive task. You're very likely to find them in variational inverse problems (especially when expressed images instead of sparse feature sets). They are very common tasks: image deconvolution, compressed sensing, dense optical flow (and more broadly dense image registration...), dense 3D reconstruction, image segmentation... ;
from my point of view, non-local image denoising (using non-local means) is very easy to understand, and it's computation-intensive enough to have been included in NVidia's Cuda examples for quite a while. It really deserves a new, better GPU implementation (such as the GPU-based kd-trees proposed a few years ago). However, it may prove challenging for an undergrad;
a useful site to see what Image Processing tasks can be done on GPUs: the Austrian (T.U. Graz) lab GPU4Vision.


Difference between Line Code and Keying The following statements are from Wikipedia page of Keying

  Keying or digital modulation applies to transmitting a digital signal
  over an analogue passband channel.

  When a digital signal is to be
  transmitted over an analogue baseband channel, the modulation
  technique is termed line coding.

To me, both of the lines seem same. The only difference I see are the terms passband channel and baseband channel.

I see both of these terms used interchangeably on the internet because we only care for the message that we get from passband channel which is obtained after filtering baseband channel.

So, what's the difference between Line Code and Keying?

I think the Wikipedia definition is limited, and anybody using the terms keying and line coding interchangably is wrong.

A line code is an encoding (or mapping) from bits to pulse amplitudes (voltages) whose main purposes are spectral shaping and receiver synchronization.

These are some of the things you can acomplish with line coding:

Remove DC. Since long transmission lines do not transfer DC very well, and since DC uses power but transmits no information, many systems use line codes that remove it. The simplest example is AMI (alternate mark inversion). In AMI, a bit 0 is mapped to 0 V, and a bit 1 is mapped to +1 V or -1 V in an alternating fashion.
Spread the signal energy across the entire bandwidth evenly. This is achieved with a scrambler.
Avoid long runs of consecutive '1's or '0's. This helps the receiver synchronize to the received signal; sync algorithmis typically rely on signal transitions. One example is the 8b10b code. This line code substitutes runs of  more than 5 symbols by sequences with transitions.

Note that the last two are different ways of achieving similar results.

It is true DC removal is more commonly used in baseband wireline systems. However, line codes, in one shape or another, are used in all communication systems.

Contrast line encoding with other types of codes: source coding compresses the information before transmission, and error control codes add extra information in order to correct and/or detect channel errors. A typical communications system will include all three types of codes.

Regarding keying: keying is the process of taking the output of the line coder and mapping it to a carrier. In this sense, keying does imply a passband channel.


Transfer function and difference equations: why does $H(z)$ numerator polynomial not correspond to $Y(z)$? For a discrete time LTI system, I understand that from a difference equation description of the system in the form
$$
\sum\limits_{k=0}^N{a_k y[n-k]}=\sum\limits_{k=0}^M{b_k x[n-k]}
$$
I can determine the transfer function of the system in the form
$$
H(z)=\frac{b_Mz^M+b_{M-1}z^{M-1}+\ldots+b_1z+b_0}{a_Nz^N+a_{N-1}z^{N-1}+\ldots+a_1z+a_0}
$$
Given a rationanal transfer function, I can easily go back to the difference equation. However, since the transfer function can also be expressed as $\frac{Y(z)}{X(z)}$,
$$
H(z)=\frac{b_Mz^m+b_{M-1}z^{M-1}+...+b_1z+b_0}{a_Nz^N+a_{N-1}z^{N-1}+...+a_1z+a_0}=\frac{Y(z)}{X(z)}
$$

Doesn't this imply that the polynomial in the numerator corresponds to the output and the denominator corresponds to the input?

But we know that the numerator polynomial derived from the input side and the denominator derived from the output side of the difference equation.

I'm obviously missing something in my understanding, but I don't see how to resolve this conundrum.  What am I missing?

It is right that $H(z)=\frac{Y(z)}{X(z)}$ and yes, the output is in the numerator. But the numerator of one side is multiplied by the denomarator of the other side to get an equality:

$$H(z)=\frac{b_Mz^M+\cdots+b_0}{a_Nz^N+\cdots+a_0}=\frac{Y(z)}{X(z)}\tag{0}$$
$$\hspace{-1.5cm}\Rightarrow\left(b_Mz^M+b_{M-1}z^{M-1}+\cdots+b_0 \right)X(z)=\left(a_Nz^N+a_{N-1}z^{N-1}+\cdots+a_0\right)Y(z)\tag{1}$$
Taking the inverse $z$ transform from $(1)$ yields the differnce equation of the system:
$$b_Mx[n+M]+b_{M-1}x[n+M-1]+\cdots+b_0=a_Ny[n+N]+\cdots+a_0\tag{2}$$
So indeed, I explained it in the reverse order. Everything begins from the difference equation of a system in $(2)$, then the $z$-transom is applied on it:
$$\hspace{-1cm}b_Mz^MX(z)+b_{M-1}z^{M-1}X(z)+\cdots+b_0 X(z)=a_Nz^NY(z)+a_{N-1}z^{N-1}Y(z)+\cdots+a_0Y(z)\tag{3}$$ and then factoring $X(z)$ on the left and $Y(z)$ on the RHS, we will get to $(1)$. From that, $(0)$ is resulted.

If we [wrongly] had
$$\frac{b_Mz^M+\cdots+b_0}{a_Nz^N+\cdots+a_0}=\frac{X(z)}{Y(z)}\tag{4}$$ then after multiplication we would  get
$$\left(b_Mz^M+b_{M-1}z^{M-1}+\cdots+b_0 \right)Y(z)=\left(a_Nz^N+a_{N-1}z^{N-1}+\cdots+a_0\right)X(z)\tag{5}$$
$$\hspace{-1.5cm}b_Mz^MY(z)+b_{M-1}z^{M-1}Y(z)+\cdots+b_0 Y(z)=a_Nz^NX(z)+a_{N-1}z^{N-1}X(z)+\cdots+a_0X(z)\tag{6}$$
and taking the inverse $z$-transform we would have
 $$\hspace{-2cm}b_My[n+M]+b_{M-1}y[n+M-1]+\cdots+b_0=a_Nx[n+N]+\cdots+a_0\tag{7}$$
which is obviously not the correct difference equation (i.e. $(2)$ that we started from).


Use Butterworth and Chebychev filters I need to calculate frequency response, phase response and apply to signals the Butterworth, Chebychev1 and Chebychev2 band-pass filters.

I'm developing in C++ with Qt, and I'm looking for algorithms or a (free) library to do that. I have my plotting solution, I just need the formulas.

The input values I have are the time/amplitude signal, the order of the filter and the low and high cut-off frequencies.

You might try Vincent Falco's Collection of Useful C++ Classes for Digital Signal Processing.

The StackOverflow Question A C++ library for IIR filter offers a few more suggestions.

Finally: the hard part for Butterworth/Chebychev filters is really the design, not the implementation.  You might consider doing the design in Matlab (or Gnu Octave), or an online site like: http://www-users.cs.york.ac.uk/~fisher/mkfilter/.  Writing the C code should be trivial once you know the coefficients for the recurrence relation.


How to segment object when it has similar histogram with background I am working in image segmentation and I would like to ask you something about my problem. I have a image such as bellow. I want to segment leopard from background. However, intensity information is not useful for my case because the bough (background) whose intensity distribution is similar with leopard's. Could you suggest to me some method to segment leopard? I attacted one example using active contour.Thank you in advance

In the first place you must admit that what you are trying to do is about impossible.

Looking at details of the image, the differences are subtle and it is very hard to tell that details belong to distinct entities.

Actually, for an algorithm there are no real "objects" to be seen but areas with different micro- and macro-texturing.

The spots on the fur are rather large and at some place you don't know if they belong to a single macro-texture or are just independent zones.

Actually, you must have previous knowledge of what a leopard is to see it. Any algorithm lacking this information will have a hard time putting the pieces together and telling the difference with the wood.

The classical approaches are by measuring texture features, such as those obtained by a bank of Gabor filters, followed by multispectral segmentation. But there is no easy roadmap, I am afraid.


What helped you in reading and implementing papers and coming up with original ideas in computer vision/graphics? I'm trying to get into academic research at the intersection of vision and graphics, and was wondering how researchers in the field went about this:

Was it useful to do a lot of fundamental math courses (Linear Algebra, Multivariable Calculus, Probability, Statistics) prior to reading papers, or did you find it more efficient to learn the math on the go?
How did you get started on reading and implementing papers? What helped and what didn't? Looking back on your experience, what did you wish you were more prepared with?

People working in computer vision/graphics may have very different backgrounds: computer programming, high-performance computing, signal processing, image analysis, optimization, machine learning, etc.

Novel contributions can be wide: from a very fast GPU/CPU algorithm to render textures (I had a reference in mind which I can't recall now) to semantic segmentation. Origins of  ideas may range from huge expertise, trials on billion of images to profound thoughts, through hardware algorithmic wisdom. Your question is far-reaching, and can dig into the depth of the psychology of scientific invention, well-described by some authors in mathematics and physics (I have Raymond Poincaré and Jacques Hadamard in mind). Here are a couple of thoughts.
The ability to (non exclusive):

be and remain curious: listen a lot, don't pretend too much, don't trust doxa (everybody knows - even yourself - that..., nobody does this), look at novelties in other fields
test and probe:hardware, data, hypothesis, models... to the limit (add too much noise, down-sample a lot, reduce learning set, cut the power),
explore unconventional ideas: transfer knowledge from one domain to the other, ask "what if I do..." (I got from Isaac Asimov that most of  science is not like Archimedes' "Eureka (I found it)",but often "Ooh, that funny"), learn XXX that strange theory, go to seminars you are not truly interesting in, or above your level
have people you trust to talk to: especially those who disagree, or dare to say "you're going nowhere", hire good colleagues, better students
take notes, talk to yourself, don't be discouraged: nice works can take decades
teach: basic student questions can drive you nuts and shake your confidence
corollary: have some solid knowledge somewhere (law of excessive learning: you can teach at one level only if you know the level above)
breathe and sleep: get away from the machine, take some fresh air, wander in woods
be patient, but remain aware of your progress: one gain confiance progressively. The first "original" ideas I got were already discovered 30 years before. I'm grateful now went I read a 1-year paper that implements what I am thinking of.
have principles: it is OK to think "This should work!" even if the evidence is shallow. Think about the "AI winter".

Those can be good incentives. On my side, I am lucky to (think that I) have a relatively broad "generic" scientific background (to connect far-fetched domains) and one or two technical skills that I dug deep into. Thus, I am not afraid of switch domains. One main risk I see is siloing: when collaborating with other, it is tempting to (and thus, try to avoid):

remain in our confidence zone,
trust others in their domain (as we would like to be treated, with respect).

One of my few scientific prides (not deep, but quite intriguing) was to work with real-time computing people (my programming skills are very limited), and discover that some "dogma" in the domain was false. My contribution was initially very small, on standard extrapolation, which I though I mastered. Some colleagues trusted me and implemented the stuff. But one colleague who did not understand the derivation asked me again, and again, and again, to understand every step that I thought was "evident". So I reduced the formulas down to the bone. With a couple of image compression tricks, we ended up with a scheme that shook the initial dogma in cyber-physical simulation. Not theoretically, but enough to feel we had "something new".


Help implementing digital filter from tansfer function I'm trying to implement a digital filter, which is given by the following transfer function:
$$
1+2V K \frac{K+c_m+2Kz^{-1}+(K-c_m)z^{-2}}{1+2Kc_m+K^2+(2K^2-2)z^{-1}+(1-2Kc_m+K^2)z^{-2}}
$$
$$
+V^2K^2\frac{1+2z^{-1}+z^{-2}}{1+2Kc_m+K^2+(2K^2-2)z^{-1}+(1-2Kc_m+K^2)z^{-2}}
$$
My interpretation of the coefficients is this:
$a0_1 =  1+2V K$
$a1_1 =  2Kc_m+K^2+(2K^2-2)$
$a2_1 = (1-2Kc_m+K^2)$
$b0_1 = K+c_m$
$b1_1 = 2K$
$b2_1 = K-c_m$
$a0_2 =  V^2K^2$
$a1_2 =  2Kc_m+K^2+(2K^2-2)$
$a2_2 = (1-2Kc_m+K^2)$
$b0_2 = V^2K^2$
$b1_2 = 2*(V^2K^2)$
$b2_2 = V^2K^2$
After dividing all the coefficients by a0, I try implementing them using a standard biquad:
$y[n] = b0*x[n] + b1 * x[n-1] + b2 * x[n-2] - a1 * y[n-1] - a2 * y[n-2]$
in parallel, with their outputs summed.
I haven't managed to get it working.
Am I on the right path? Is there anything I'm doing terribly wrong, or not doing? Thanks very much for any help!

You seem to do a few things wrong

Start with just transcribing each fraction, ignore the factors in front of the fraction and "1" in front
Make sure you sort by powers of $z$
Normalize all coefficients to $a_0$
Multiply all $b$ coefficients with the factor for the fraction.
The "1" is a third filter

Example
$$a0_1 = 1+2Kc_m+K^2 \\ a1_1 = 2K^2-2$$


How to merge 2 WAV files together? Can I take the average of the 2 audio file float values? I'm trying to do one simple thing on Google's C++ Oboe library: pick 2 recordings and merge both. That is, make as if the 2 were playing at the same time.
Oboe works with Float numbers from -1.0 to +1.0. Can I simply do like this:
Audio1:
-0.2, -0.221, -0.23, -0.24, ...

Audio2:
-0.3, -0.23, -0.12, -0.87, ...

Audio Mixed:
(-0.2+-0.3)/2,  (-0.221+-.023)/2, (-0.23 -0.12)/2, (-0.24-0.87)/2, ...

?
That is, take the average of both audios into the new audio?
Is this how 2 simple audio files are merged together?

If your goal is to simply make it sound like they were playing at the same time, you can simply just add them together element-wise. There is need to "average". The resulting waveform is the same, but scaled by a factor of 2 in your case. Technically, that number could really be anything.
If you want different properties, where maybe one sounds louder than the other, then you'll have to do something else.