Leakage in power spectral density estimation What does it mean for a power spectral density exhibiting leakes from one frequency bin to other.

I am reading a book which states the following:

"Often the time series we use as input for our algorithm will have a non-zero DC average, which
may even slowly change over time. Such a DC average will show up in the first frequency bin
(m = 0) of the resulting spectrum. If a window function is used or the average changes over
time, it will also leak into adjacent frequency bins, possibly masking low-frequency signals."

Could somebody explain what does this mean? Why does a window leak power into adjacent frequency bins?

I also read the following which has probably a similar connection and it would be good to get an explanation on this as well:

"If we simply take a stretch of length N out of a time series containing a sinusoidal signal and
perform a DFT, we will most likely find that the sinusoidal signal which we might naively expect
to result in a sharp peak in only one frequency bin, will instead show up as something ugly, The reason is that the DFT implicitly assumes that the signal is periodic, i.e. that the time series
of length N repeats itself infinitely in a cyclic manner. If the frequency of the sinusoidal input
signal is not an exact multiple of the frequency resolution fres, i.e. does not fall in the exact
center of a frequency bin, this assumption is not true, and the DFT will ‘see’ a discontinuity
between the last sample and the first sample due to the cyclic continuation. That discontinuity
spreads power all across the spectrum"

What does it mean for DFT to "see" a discontinuity as written above"?

EDIT:After accepting answer:

I ran a MATLAB simulation (code below) in which the frequency of sinusoidal is exactly a multiple of resolution of the PSD and thus falls directly at the center of a frequency bin, still there is leakage into neighboring bins (3 to the right and three to the left of this bin in which the sinusoidal lies). This spilling gets more and more resolved as I increase the window length. Which would seem to suggest that even the frequencies that fall directly at the center of frequency bins leak into the adjacent frequency bins and the amount of leakage will depend on the shape of the window and the length of the window. Is this understanding correct?

MATLAB code:

%%%%%%%%%%%%%%%%%%%% PSD Estimation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all;
N=50; %number of periods

Fs= 4096;
snr=10;

M=256; %length of FFT

Fres=Fs/M;

Mul_factor=20;

F_c=Fres*Mul_factor;  %in Hz%

Window_length=128;

Overlap=floor(0.4*Window_length);
x=0:1/Fs:N*(1/F_c);
y=sin(2*pi*x*F_c) + 1/sqrt(snr).*randn(1,length(x));

[ppx,w]=pwelch(y,Window_length,Overlap,M,Fs);

figure;
hold on;
subplot(2,1,1);
plot(x,y);
subplot(2,1,2);
plot(w,ppx);
hold off;

Could somebody explain what does this mean? Why does a window leak
  power into adjacent frequency bins?

You can think of frequency bins as imposing an artificial frequency grid which allows only certain frequencies to fit nicely within it (producing a sharp peak) and all other frequencies are smeared across adjacent bins thus producing spectral leakage.

This means that only those frequencies that are centered exactly on bin frequency will not produce leakage.

Since the bandwidth of a bin (i.e frequency resolution) is simply the
sampling rate divided by the dft length, we can conclude that only frequencies that are exact integer multiples of the frequency resolution will not produce spectral leakage.

For example, if the sampling frequency is 8000 Hz and dft length is 1024, the frequency resolution is computed as

 frequency_resolution = 8000 / 1024 = 7,8125 Hz

So, if your signal contains frequencies which are integer multiples of 7,8125, you'll get
no leakage (e.g 1*7,8125, 2*7,8125, 3*7,8125 and so on)

If, however, your signal contains a frequency that is somewhere between 2 adjacent bins, you'll get spectral leakage, which I think is intuitive because no single bin can
fit that frequency component exactly.

UPDATE (trying to answer the comments below):

It's true that windowing does introduce some additional leakage on its own. But even if you don't use windowing at all (which corresponds to a rectangular window), you'll still get spectral leakage (provided that the frequency under examination is not an integer multiple of the frequency resolution).
In fact, since using a rectangular window (no window at all) produces the worst leakage, using overlapping windows other than the rectangular helps reduce the spectral leakage.

So, why do we have spectral leakage in the first place?

Remember that DFT assumes that a signal is periodic. Since the frequency resolution is  fixed (e.g 7,8125 Hz in the example above) for every bin, which in this case corresponds to exactly 1024 samples, we can see that a single bin can hold exactly one cycle of a 7,8125 Hz frequency , 2 cycles of 15,625 Hz frequency, 3 cycles of 24,4375Hz frequency and so on. This is true for every bin, which means that there'll be no discontinuities either within a bin or between adjacent bins.

However, if we now change our frequency to something that fits between two bins, things change. This new frequency will take either more (or fewer) than 1024 samples (one or more cycles), which means that we'll get discontinuities either within a bin or between neighbouring bins and this in turn creates spectral leakage.

Hope this helps.


Please clarify the mandatory condition on Vestigial Sideband Modulation filter? In VSB filter we have the Transfer function condition as $H(F-F_c)+H(F+F_c)=1$
So can we accept that for the range of $F$ in $(F-df,F_c)$ we have $H(F)$ and for $(F_c,F+dF)$ we have $A-H(-F)$  where $A=H(F_c)$.
I am learning Latex so need some time:(

For $(F_c,F+dF)$, it needs to be $2A-H(-F)$, when $H(F_c)=A$. See the picture below taken from (http://163.152.6.72/contents/vod/21/20121129093403/20121129093403.pdf). Here $A$ corrsponds to $1/2$ as per below figure, that is the value at $H(F_c)$. In order for the sum to become $1$, $A$ needs to $1/2$. Also, $\epsilon$ is value at any frequency in the range $(F-dF,F_c)$. So the value in $(F_C,F+dF)$ needs to be $2A-H(-F)$ so that sum $H(F) + 2A-H(-F) = 1$, as $H(F)=H(-F)$ along $F=0$.


How does one compute the band width of a signal in the time domain $x(t)$? Sorry if this is a super easy question, I am really not a signals guy (more of a statistician, computer scientist). But I was trying to reconstruct/learn a sinusoidal function $x(t)$ with linear regression $f(x) = \langle w , \phi(x) \rangle $ (I'm adding as many polynomial features as I may need) over a fixed interval $[-1,1]$. Then I realized that there is this thing called the Nyquist-Shannon sampling theorem that might help figure out if I have enough samples. The reason that I care about it is because the performance of my regression model is really bad on the test set and I wondering if I am just not getting enough samples from my ground truth signal (which for the moment I control synthetically). I know now that I need my sample frequency $f_s$ to be larger than twice the band width $B$ of my signal:

$$ f_s > 2B$$

to my understanding bandwidth is just the range of the frequencies present in the signal $x(t)$ (I even asked to clarify that to me just in case here: What is the definition of the bandwidth of a signal?). Thus if that is correct then band width is:

$$ B = f_{max} - f_{min} $$

is there a general way to get $f_{max}$,$f_{min}$ from a signal in the time domain?

If I understand this correctly then I just need to get the frequencies that define $x(t)$? Right?

If that is correct then the only way I know is if we know the analytic form of the $x(t)$, which I do for my example, then I just read of the frequencies from the Fourier series $ f(x) = \frac{1}{2}a_0 + \sum^{\infty}_{n=1} a_n cos(nx) + \sum^{\infty}_{n=1} b_n sin(nx) $. Say $n_{max},n_{min}$ are the largest values modifying the inside of sin or cos terms then we read them off and do:

$$ B = \frac{n_{max} }{2 \pi} - \frac{n_{min}}{2 \pi} $$

right? What happens if there is only 1 single term like it only a sin like $x(t) = sin(n t)$, how do we get the smaller frequency if there is no other frequency?

Now that I understand much better what the term frequency means I assume that the answer will point me that if I don't know the analytic form of my signal (as probably happens most in practice) we need to resort to some type of transform used. Probably the Fourier transform? But a quick google search to the frequency domain yielded a list of methods:

Fourier series – repetitive signals, oscillating systems
Fourier transform – nonrepetitive signals, transients
Laplace transform – electronic circuits and control systems
Z transform – discrete-time signals, digital signal processing
Wavelet transform - image analysis, data compression

I don't know when to apply each but it seems the wikipedia article links to more articles for each one...

Is this correct or am I totally off? (which might be the case since I'm not very familiar with this field). I am assuming there are probably engineering details I am totally missing like flow pass filters or something else that might be important for real signals in practice. Though its just a random guess based on the definition of low-pass filter (filters high frequencies out). Though not sure.

Further reading What is the definition of the bandwidth of a signal? makes me suspect that the actual correct answer is as follows:

Consider a time domain signal $x(t)$.
Consider its frequency description $X(\omega)$

the the support of $X(\omega)$ is the bandwidth. In other words, "the summations" of the points in the domain (to be more precise, the integral of the domain). So if the domain is the whole interval $[-W,W]$ the the bandwidth would be:

$$ B = \int^{W}_{-W}d \omega = 2W $$

so in general is it:

$$ B = \int_{w \in supp(X(\omega))} d\omega $$

would be my guess.

As far as I know there is no way to get $f_{max}$ and $f_{min}$ directly from time domain (zero crossing rate will not help if the signal has a non-zero DC offset for example) you need to look at the frequency domain.

If there is just one term, the smaller frequency is just $0$

Given a time signal and its sampling rate ($Fs$), you can take the discrete Fourier Transform to study its frequency spectrum. Unfortunately you cannot use this to decide your sampling rate if your initial sampling rate ($Fs$) is lesser than the 2$*$bandwidth (2$*$[highest frequency] if you are looking for perfect reconstruction) as you have already lost the data during initial sampling(look up aliasing for further reading).

Bottom Line: You need to have an estimate in what frequency range your signal will lie so that you could sample it accordingly.


Precise Centre frequency of an All-pole digital filter I'm not an engineer and have essentially taught myself all I know - this present problem is giving me some problems.

I have an all-pole filter (a gammatone filter - impulse response is essentially a damped cosine) implemented as a cascade of identical 2nd-order pole pair sections. The poles are located in the z-domain at $\exp(-\beta + i\theta)$: $\beta$ is the bandwidth in radians/sample and $\theta$ is the 'ringing' frequency in radians/sample.

Given a value for $\theta$, this pole pair does not produce a filter with peak gain at the frequency $\theta$ - I'm trying to find out exactly what the offset is (i.e. at what frequency a filter with a given $\theta$ actually has peak magnitude), so that I can produce filters that peak at precisely a given arbitrary frequency.

I've found an equation for calculating what this offset is given poles located in the s-domain at $-\beta + i\theta$:

$$\mathrm{f_{c,actual}} = \theta  \sqrt{1 - \dfrac{\beta^2}{\theta^2}}$$

This formula provides results consistent with observations up to a certain frequency, though I believe it needs to be modified to produce results that work for arbitrary discrete spectrums.. I've tried a number of things and just can't quite figure it out.

Any help is greatly appreciated.

For a filter consisting of a complex conjugate pair of poles, the $z$-domain transfer function is:
$$H(z) = \frac{a}{\left(1-r(\cos{\theta}-i\sin{\theta})z^{-1}\right)\left(1-r(\cos{\theta}+i\sin{\theta})z^{-1}\right)}\\
= \frac{a}{1 - 2r\cos(\theta)z^{-1} + r^2z^{-2}},$$

where $a$ is a constant by which the magnitude frequency response can be normalized, $r$ is the pole radius, and $\theta$ is the pole angle. The peak of the magnitude frequency response will not be exactly at frequency $\theta$, but at a nearby frequency $\omega_0$ at which the derivative

$$\frac{|a|2r\sin(\omega)\left(2r\cos\omega - (r^2 + 1)\cos\theta\right)}{\left(4r^2\cos^2\theta - 4r(r^2 + 1)\cos(\theta)\cos(\omega) + 4r^2\cos^2\omega + r^4 - 2r^2 + 1\right)^{3/2}}$$

of the magnitude frequency response $|H(e^{i\omega})|$ with respect to $\omega$ equals 0. That frequency is:

$$ω_0 = \operatorname{acos}\left(\frac{(r^2 + 1)\cos\theta}{2r}\right)\text{ if }-1 < (r^2 + 1)\cos(\theta)/(2r) < 1$$

Figure 1. The condition for that the magnitude frequency response peaks elsewhere than at $w_0=0$ or $w_0=\pi$ is that poles stay outside these droplets.

In the following it is assumed that the condition is satisfied as otherwise the peak would be trivially at frequency $0$ (if $\theta < \frac{\pi}{2}$) or $\pi$, whichever is closer to the pole. The normalized magnitude frequency response with $a$ chosen so that $|H(e^{i\omega_0})| = 1$ reads:

$$|H(e^{i\omega})| = \sqrt{\frac{(r^4 - 2r^2 + 1)\sin^2\theta}{4r^2\cos^2\theta - 4r(r^2 + 1)\cos(\theta)\cos(\omega) + 4r^2\cos^2\omega + r^4 - 2r^2 + 1}}.$$

Given $\omega_0$ and $r$:

$$\theta = \operatorname{acos}\left(\frac{2r\cos\omega_0}{r^2 + 1}\right)\tag{1}$$

which is something you can use directly to keep the peak of your all-pole bandpass filter where you want it:

$\omega_0$ pole loci">
Figure 2. Pole loci for a selection of constant $\omega_0$ traced over $0 < r < 1$.

Applying the substitution given by equation (1) gives:

$$H(z) = \frac{a}{1 - \frac{4r^2\cos(\omega_0)}{r^2 + 1}z^{-1} + r^2z^{-2}}$$

Just like $\theta$ does not exactly give the peak frequency, $r$ does not give the exact bandwidth. The exact -3 dB to -3 dB bandwidth $B_\text{-3 dB}$:

Figure 3. -3 dB bandwidth of the filter, with $\omega_0 = \frac{\pi}{4}$ and $r=0.9$.

is given by:

$$B_\text{-3 dB} = \operatorname{acos}\left(\frac{(r + 1/r)\cos\theta + (r - 1/r)\sin\theta}{2}\right)\\- \operatorname{acos}\left(\frac{(r + 1/r)\cos\theta + (1/r - r)\sin\theta}{2}\right)$$

This was found by solving the values of $\omega$ for which $|H(e^{i\omega})| = \sqrt{2}/2$. The first $\operatorname{acos}$ gives the solution found to the right from the peak and the latter $\operatorname{acos}$ the one to the left. More generally, for a bandwidth from magnitude frequency response value $g$ to $g$:

$$B_g = \operatorname{acos}\left(\frac{(r + 1/r)\cos\theta + \sqrt{1/g^2 - 1}(r - 1/r)\sin\theta}{2}\right)\\-\operatorname{acos}\left(\frac{(r + 1/r)\cos\theta + \sqrt{1/g^2 - 1}(1/r - r)\sin\theta}{2}\right)$$

This is needed when designing a composite filter, for example a gammatone filter, that is a cascade of identical two-pole filters. Usually we want to design the filter to have a given bandwidth, but it is very difficult to solve $r$ from the above. So we resort to approximation. Using the value $D_2$ of the second derivative of the magnitude frequency response at $\omega_0$:

$$D_2 = \frac{(r^2 + 1)^2\cos^2\theta - 4r^2}{(r^2 - 1)^2(1 - \cos^2\theta)},\tag{2}$$

we can approximate the magnitude frequency response with that of a complex one-pole filter that has at $\omega_0$ the magnitude frequency response and its first two derivatives matched to those of the two-pole filter:

Figure 4. One-pole approximation with pole angle $\omega_0$.

The magnitude frequency response of the one-pole filter is:

$$|\hat H(e^{i\omega})| = \frac{\sqrt{2}\left(\sqrt{1 - 4D_2} - 1\right)}{2\sqrt{\left(2D_2\cos(\omega - \omega_0) - 2D_2 + 1\right)\left(1 - \sqrt{1 - 4D_2} - 2D_2\right)}}.$$

It is symmetrical around $\omega_0$ and has $g$ to $g$ bandwidth:

$$\hat B_g = 2\operatorname{acos}\left(\frac{1/g^2 - 1}{2D_2} + 1\right),$$

or inversely and more usefully:

$$D_2 = \frac{1 - g}{2g^2\left(\cos(\hat B_g/2) - 1\right)}.$$

Equations $(1)$ and $(2)$ when combined give:

$$D_2 = \frac{2r^2(r^2 + 1)^2\left(1 - \cos(2\omega_0)\right)}{(r^2 - 1)^2(2r^2\cos(2\omega_0) - r^4 - 1)}.$$

Solving that for $r$ gives a cumbersome equation with $\omega_0$ and $D_2$ as arguments that can be used to set the approximate bandwidth of the two-pole filter (c language compatible plaintext with w_0 = $\omega_0$ and d_2 = $D_2$):

r=sqrt(-sqrt(2)*sqrt(sin(w_0))*sqrt((-sin(w_0)*(pow(d_2+1,2)*pow(cos(w_0),2)
+d_2-1)*sqrt(-pow(d_2+1,2)*pow(cos(w_0),2)+pow(d_2,2)-2*d_2+1)+pow(d_2+1,3)
*pow(cos(w_0),4)-(d_2+1)*(pow(d_2,2)-d_2+2)*pow(cos(w_0),2)+pow(d_2,2)-2*d_2
+1)/sqrt(-pow(d_2+1,2)*pow(cos(w_0),2)+pow(d_2,2)-2*d_2+1))+sin(w_0)*sqrt(-
pow(d_2+1,2)*pow(cos(w_0),2)+pow(d_2,2)-2*d_2+1)-(d_2+1)*pow(cos(w_0),2)+1)
/sqrt(-d_2)

or with alternative symbols:

r = √(- √2·√(SIN(ω_0))·√((- SIN(ω_0)·((d_2 + 1)^2·COS(ω_0)^2 + d_2 - 1)·√(-
(d_2 + 1)^2·COS(ω_0)^2 + d_2^2 - 2·d_2 + 1) + (d_2 + 1)^3·COS(ω_0)^4 - (d_2
+ 1)·(d_2^2 - d_2 + 2)·COS(ω_0)^2 + d_2^2 - 2·d_2 + 1)/√(- (d_2 + 1)^2
·COS(ω_0)^2 + d_2^2 - 2·d_2 + 1)) + SIN(ω_0)·√(- (d_2 + 1)^2·COS(ω_0)^2 +
d_2^2 - 2·d_2 + 1) - (d_2 + 1)·COS(ω_0)^2 + 1)/√(-d_2)

When $\hat B_\text{-3 dB} = \hat B_{\sqrt{2}/2}$ is set proportional to $\omega_0$ by a constant $c$ stepped from 0.1 to 1.0, the upper half-plane pole resides on tracks like these:

Figure 5. Pole loci for bandwidth proportional to peak frequency with proportionality constant $c$.


Removing noise from audio using Fourier transform in Matlab I'm trying to remove noise from an audio file. This audio file contains speech as well as constant pink noise. I know that I have to use the Fourier transform to convert to the frequency domain then use a filter to filter out the frequencies of the pink noise, but I really don't know where to start with that. It also seems that the noise covers many more frequencies than the speech file (I have the speech and noise files separately) by looking at the plots in the frequency domain.

Doesn't that mean that if I filter out the frequencies of the noise it'll remove all the speech as well?

The plot of the speech audio in the frequency domain displays that the frequencies are mostly between -450 and 450.

The plot of the noise alone in the frequency domain displays that the frequencies are between -2000 and 2000. This clearly covers a much wider range of frequencies than the speech audio.

Does this mean I can't remove it from the noisy speech audio without removing the speech too? Or can I remove the frequencies outside of -450 and 450 and that would remove much of the noise? If so, how do I compute the range of frequencies instead of looking at plots myself to see the range?

From my experience, your best bet is spectral subtraction. You can find a lot of literature about it on the web. Here just a quick overview of what you need to do. You take overlapping windowed blocks of your time domain signal and transform them to the frequency domain using an FFT. Then for each frequency bin you need to estimate the signal-to-noise ratio. There are several simple noise tracking algorithms that perform well if the noise is relatively stationary. Based on the estimated SNR per frequency bin you mulitply each bin with a gain constant between 0 (terrible SNR) and 1 (no noise). Note that you only deal with magnitudes, the phase of the FFT is left unchanged. After modifying the FFT bins with the gain factors, you combine the processed bins back to the time domain using the phase of the original noisy signal.

Here is a classic paper on this subject. It goes a bit further than pure spectral subtraction but the overall system is the same. Don't worry too much about the math involved to derive the optimal gain function.


Why are Gaussian filters used as low pass filters in image processing? In 1d signal processing, many types of low pass filters are used. Gaussian filters are almost never used, though.

Why are they so popular in image processing applications? Are these filters a result of optimizing any criterion or are just ad hoc solution since image 'bandwidth' is usually not well defined.

Image processing applications are different from say audio processing applications, because many of them are tuned for the eye. Gaussian masks nearly perfectly simulate optical blur (see also point spread functions). In any image processing application oriented at artistic production, Gaussian filters are used for blurring by default.

Another important quantitative property of Gaussian filters is that they're everywhere non-negative. This is important because most 1D signals vary about 0 ($x \in  \mathbb{R}$) and can have either positive or negative values. Images are different in the sense that all values of an image are non-negative ($x \in  \mathbb{R}^+$). Convolution with a Gaussian kernel (filter) guarantees a non-negative result, so such function maps non-negative values to other non-negative values ($f: \mathbb{R}^+ \rightarrow \mathbb{R}^+$). The result is therefore always another valid image.

In general, frequency rejection in Image processing in not as crucial as in 1D signals. For example, in modulation schemes your filters need to be very precise to reject other channels transmitted on different carrier frequencies, and so on. I can't think of anything just as constraining for image processing problems.


FIR filter : $h(n)$ is palindrome? I study condition for fix if filter is a linear phase,but it's not clear in my mind!

I have this $h(n)$:
$$h(n) = \begin{cases} \left(\frac{1}{2}\right)^{n} & 0<n<N-1 \\[2ex] 0 & \text{elsewhere} \end{cases}$$
I got the $H(z)$:
$$
H(z) = \frac {\left(\frac{1}{2}\right)^{N}{z}^{-N}-1}{\frac{1}{2}{z}^{-1}-1}
$$
Now, for the phase linear what I do ? I think to verify palindrome property for this $h(n)$, but how ?? or I can finally say that to the form that has definitely not linear ?

There are four types of linear phase FIR filters, as described in this answer. They have to satisfy the following symmetry conditions:

$$h[n]=h[N-1-n]\quad \text{or}\quad h[n]=-h[N-1-n]\tag{1}$$

(assuming the impulse response $h[n]$ is zero for $n<0$ and $n\ge N$). Since your filter satisfies neither of the two conditions given by $(1)$ it can't have a linear phase response.


Square of Frequency characteristic of the system What does $|H(j\omega)|^2$ in $20 \log_{10}|H(j\omega)|$ mean?

Is it some ratio of energy or power? And why? How to derive it?

I'm sorry for curtness of my questin.

I suppose you're talking about the frequency response of a continuous-time linear time-invariant system. It is defined as the Fourier transform of the system's impulse response $h(t)$:

$$H(j\omega)=\int_{-\infty}^{\infty}h(t)e^{-j\omega t}dt$$

Its practical relevance is that it shows the frequency dependence of the system's input-output relation. If $X(j\omega)$ is the Fourier transform of the input signal, and $Y(j\omega)$ is the Fourier transform of the output signal, then
$$H(j\omega)=\frac{Y(j\omega)}{X(j\omega)}$$

For a sinusoidal input signal $x(t)=\sin (\omega_0 t)$ the output signal $y(t)$ is given by

$$y(t)=|H(j\omega_0)|\sin (\omega_0 t + \text{arg}\{H(j\omega_0\})$$

where $\text{arg}\{H(j\omega_0\})$ is the phase of the complex function $H(j\omega)$ at frequency $\omega_0$. So you can see that $|H(j\omega)|$ represents the amplification of the system at frequency $\omega$. From above relations it is also obvious that the output signal cannot contain any frequencies that were not present in the input signal. $H(j\omega)$ can only reduce or amplify frequency components of the input signal. This is how (linear time-invariant) filtering works.


Is there any way to measure of Gaussian-ness? I have some sampled data that has $1/f$ noise in it, with departures from the mean. These are long term departures. I could use something like a median filter but the window length would be longer than I would like.

Is there any way to measure its departure from a Gaussian distribution?
Or measure how 'Gaussian-like' a time series statstical sample is?

One way to measure the similarity between two distributions is the Kullbeck-Leibler divergence. Granted, it is not a norm because it is not symmetric, but it is a way of quantifying the distance between two pdfs.

EDIT: There is a topic called "Normality testing" or "Gaussianity testing" that is devoted to that very task. Perhaps some tools from that field will help.


how important is the pre-emphasis stage in MFCC? I wrote a simple MFCC routine, using all the steps usually suggested in the literature with the exception of the pre-emphasis stage. How important is this stage?
 Also, some authors suggest that using rectangular windows instead of triangular doesn't degrade MFCC performance (as long as the windows are overlapping) so any ideas regarding this particular subject are welcome.

Pre-emphasis is a way of compensating for the rapid decaying spectrum of speech. The experiment is worth trying on real data - you will find that the DCT basis is better at extracting a set of decorrelated coefficients when the spectrum has been whitened by the pre-emphasis filter.

This justification does not really apply in the case of music where the spectrum is more slowly decaying - it is not uncommon in MIR papers to see MFCC computed without the pre-emphasis filter. In the case of music, an argument for pre-emphasis is that it provides a very crude approximation of the outer ear filtering. Some feature extraction modules are indeed using MFCC-like features with a proper outer ear filter (response given by Terhard's formula) in place of pre-emphasis - For example, it looks like this is what is used by the Echo Nest analyzer.


Why is $x$ not considered a primitive polynomial while being considered an irreducible polynomial? In "Shu Lin, Daniel J. Costello-Error Control Coding (2nd Edition), Prentice Hall 2004" it is given that in $GF(2)$, if $f(x)$ is an irreducible polynomial of degree $m$ it divides $x^n+1$ where $n = 2^m -1$. And as per mathematica documentation, $x$ is an irreducible polynomial, here is the conflict.  $x$ does not divide $x+1$.  Also $x$ is declared non-primitive!  Just want clarification.

From the Wikipedia page irreducible polynomials: "a non-constant polynomial $f$ in $F[x]$ is said to be irreducible over [the field] $F$ if it is not the product of two polynomials of positive degree."

From this definition, it should be clear that $x$ cannot be viewed as the product of two polynomials of positive degree. Take one polynomial to be $x$ then there is no other polynomial (with degree greater than zero). Irreducible polynomials act like prime numbers in the field in that other polynomials can be factored into (powers of) the irreducible polynomials.

A primitive polynomial has a different definition. "In field theory, a branch of mathematics, a primitive polynomial is the minimal polynomial of a primitive element of the finite extension field $GF(p^m)$." Notice that irreducible polynomials can be defined over any field, but primitive polynomials only apply to the extension of a finite field. These polynomials are minimal polynomials that generate all of the elements of the extension field.

From this definition, we require that $x$ generates all of the (nonzero) elements of the field. It does in this case. We also require that it is a minimal polynomial which means it has a root which is a primitive element (generates all of the nonzero elements of the extension field). It fails this test as the root of $x$ is $0$ and $0$ is not a primitive element of any field since $0^n = 0$.


Is there a general relation between mean power and period I posted this question at Math Stack Exchange a few weeks ago without any luck.  It occurred to me I might have more luck by posting it here.  The basic question is whether there exists a known relationship between the period of a signal (i.e. the inverse of the frequency) and its corresponding power, such as that depicted in the graph in the linked post.

There is no universal relationship between power and period length. These are two unrelated, yet completely independent, concepts (see amplitude vs. frequency).

But what your plot might be referring to is a digital CMOS circuit. Such a circuit consumes a certain amount of energy per clock cycle. Hence for this kind of device, power consumption is (approximately) proportional to clock frequency.

So it depends on the origin of the data. There is no universal relation, but there might be one for specific situations.


In image compression using transforms, how to deal with the transformed coefficients as they are not integers? I am new to the field of image compression. While going through various texts, I read about how transforming the image to another domain using, for example, the wavelet transform, or the DCT, makes it suitable for lossy compression.

I have been trying to play around with wavelet transforms on a sample bmp in MATLAB. Using the MATLAB library function, I created transform coefficients using Daubechies wavelets.

Now a given pixel when read from bmp is in uint8 format, i.e., it uses up 1 byte of memory. This is because in most image formats the range of values a pixel can take is from 0 to 255. After wavelet transformation however, one gets double precision floating point for each element of the transformed matrix. This would typically consume about 4 bytes for each element.

My goal was to try out some very basic compression using quantization and entropy coding of quantized elements. Quantizing the low magnitude values of the transformed elements to zero does give good compression. But trouble is I have to still use double precision floating point. The moment I convert all the wavelet coefficients to some integer format, the reconstructed image becomes very poor in quality as compared to the original.

In standard compression techniques, how are the transformed coefficients handled, given that they are not integer in format? Each transformed coefficient, due to its double precision nature, takes up more bytes than the original. Can you point out resources where this issue has been explained?

Oftentimes some rounding occurs in storing the coefficients. This is why many image compression algorithms are lossy, i.e. they lose information when converting the floating point coefficients to integer format. The process of rounding is called quantization. See this wikipedia article for an example.

http://en.wikipedia.org/wiki/JPEG#Quantization


Base line detection of text lines I came across various works in the literature of document analysis addressing the problem of detecting base lines in multi-oriented texts. But no amount of search could give me satisfactory answer to the following quest.

Is there an algorithm to detect the base line which neither use a histogram based method or any kind of threshold parameter?

I don't think you can completely get rid of histogram or threshold-based binarization, since the former is to achieve line segmentation, while the latter is to extract the letters. The Radon horizontal projection is used for line segmentation, and the center line can be used to approximate the baseline of each segment. Yet this is somehow equivalent to the horizontal histogram.

There are some improvement based on the basic histogram. In Zheng et al's paper  A model based line detection in documents, they applied a directional single connected chain strategy that considers the global textline orientation and merges the non-overlapped small segments into one line. The algorithm provides more accurate base line detection than the pure histogram-based line segmentation.

Besides, in Arabic script where short characters and large diacritics exist, the Tarik et al's paper A Novel Baseline Detection Method of Handwritten Arabic-Script Documents Based on Sub-Words that measures the vertical distance between the semi-overlapped subwords baselines.

Yet in those studies in the attempt to increase the accuracy, histogram and/or threshold were still a necessity.


Channel impulse response in freqency domain I have a doubt about DFT matrix and channel impulse response.

Let $h$ is a channel impulse response , $F$ is a DFT matrix. IF I take a product. I got:
$$ F* h=\tilde{h}$$

$\tilde{h}$ is a frequency response of the channel, isn't? or channel impulse response in frequnce domain

You are right. If you have MATLAB, you can do a little experiment to double check yourself:

h = randn(100, 1);     % make some random impulse response
F = dftmtx(length(h)); % make the DFT matrix (100-by-100)
abs(fft(h) - F*h)      % this number should be small


Component of one signal along the another Whenever we need to find component of one signal along another signal we can do so by the help of dot product. I want to as if this method is valid for all kind of signals. Are there any restrictions upon either of the signals?

Generally speaking, about the only restriction that the signals need satisfy is that they have finite energy or power.  Since you mention vectors by which is presumably meant $N$-tuples of real or complex numbers and dot products, the restriction is the trivial one that each component of the $N$-tuple has finite value.

Practically, there are several issues that you need to think about depending
on what you are trying to do.  Obviously
$$\langle \mathbf x, \mathbf y \rangle = \sum_{n=1}^N x[n](y[n])^*$$
is the dot product but how accurate a measure it is of the inner product
$$\langle x(t), y(t) \rangle = \int_{-\infty}^{\infty} x(t)(y(t))^* \mathrm dt$$
of the continuous-time signals $x(t)$ and $y(t)$ from which these samples
came depends on whether the signals were filtered before sampling, the
sampling rate, whether the samples span the entire time interval where
the signals are nonzero, whether windowing was applied to the samples, etc.
You might be interested in the cross-correlation function $R_{x,y}(\cdot)$
instead of just the inner product (which is $R_{x,y}(0)$) in which case
FFTs might be used and you have to make sure that the signals are properly
zero-padded so that you get the right answers. Finally, we began with
the "trivial" restriction that each component is finite which always holds
in practice, but if the components are large in magnitude, you need to worry
about  the possibility of overflow in the inner product computation.

So, as you can see,
in theory, theory and practice are the same, but in practice, they are
different.


Proof that DFT does not require more than N points I'm trying to show how the discrete Fourier transform (DFT) arises from the equation for the continuous-time Fourier Transform. I've run into an interesting caveat which I can't seem to find an explanation to anywhere. My 'derivation' goes something like this:

Lets consider the continuous time Fourier transform equation
$$
F(\omega) = \int_{-\infty}^{+\infty} \, f(t) \, e^{-2\pi i \omega t} \,dt
$$
If $f(t)$ was to be discretised at $N$ evenly spaced sampling points with the sampling interval of $dt$, then it could be thought of as a continuous function of time with Dirac deltas at each sampling point. The integral of each of the deltas would be equal to the value of $f(t)$ at that point. Lets denote $f(t)$ at $k$-th sampling point (0-indexed) as $f_k$. Now, the integral is non-zero only at those delta functions so the entire above equation can be written as a sum
$$
F(\omega) =  dt \times \sum_{k=0}^{N-1} \, f_k\, e^{-2\pi i \omega (k dt)}
$$
Even though the sum is finite, this is a continuous function of $\omega$. $F(\omega)$ turns out to be periodic in $\omega$ and therefore just a single period of it is needed to be kept to retain all of the information about the  sampled version of $f(t)$ it came from.

All of this is (I believe) fine. Here comes the tricky bit. I want to then say that this single period of $F(\omega)$ can be sampled down to $N$ evenly spaced points which would in fact give the equation for the DFT. I say that: "It is one of the implications of the Shannon information theory (...) that for a general case of $N$ arbitrary bits of information, there does not exist a smaller set of bits which can be used to represent all of the information of the original $N$ bits. This places a lower bound on how many samples of $F(\omega)$ one can take."

Question 1: Is this a correct statement about implications of the information theory?

Question 2: I know that N bits of $F(\omega)$ are sufficient, but how do I prove it? To put it differently: How can I one show that for any arbitrary signal I will not loose any information about $f(t)$ by sampling $F(\omega)$ at only $N$ points.

Note that f() can only be sampled down to N points per period if f() is bandlimited to include only spectrum below t_period/(2*N) in frequency.  And a signal exactly periodic in N, with no frequency content at or above t_period/(2*N) in frequency can be exactly represented with N complex exponential basis vectors (Fourier).  And the inverse.

The number of degrees of freedom in N complex basis vectors or N complex DFT results matches the number of degrees of freedom in N complex time samples.


Estimating the Variance of Noise Using Median Function Applies on the Gradient of the Image By investigating the implemented code of an accepted paper, I encounter with following relation to estimate the variance of noise in an image:

$$ \sigma = \frac{median (|\nabla^hx-median(\nabla^hx)|)}{2*.6745},\tag{1}$$
where $\nabla^h$ compute horizontal derivative of $x$.

I know the variance of some samples can be estimated by $\frac{1}{n}\Sigma_{i=1}^n (x_i-\mu)^2$, but I cannot understand $(1)$ as well as the magic number $.6745$.
Any help will be appreciated.

This method is called median absolute deviation (MAD). Unlike the version with mean, it is a robust estimator of the variance. For a univariate set $X$, MAD is defined as:

$MAD = median_{i} (| X_i-median_j(X_j) |)$.

This is typically thought as the median of the absolute values of the deviations from the median of the data. The variance is then robustly estimated as:

$\sigma= \kappa MAD(X)$, where $\kappa=1.4826$ is a scale factor. Note that your $0.6745=\frac{1}{1.4826}$ corresponds to such factor.

In your example, MAD of the gradients are utilized to estimate variance of the noise. It is a pretty standard procedure. For more on the topic, please refer here.

In MATLAB, you could simply compute this via the function mad.


Information Theory - units of channel capacity In a first course in Information Theory, when the operational interpretation of channel capacity is introduced, it is said to be the highest data rate (in bits/channel-use) of reliable communication. While reading a few papers, I came across channel capacity being expressed in units of bits/s/Hz. So I was thinking about the connection between the two units and came up with the following explanation. Please let me know if this is wrong.

For a bandlimited channel (bandwidth = $W$ Hz), you can transmit at $2W$ symbols/sec by the Nyquist sampling theorem. So the rate "per bandwidth" (spectral efficiency) can be written as 2 symbols/sec/Hz. If each symbol is 1 bit, then you are transmitting 1 bit in each of the samples. So is 1 bit/channel-use equivalent to 2 bits/sec/Hz?

What is one "channel-use"?

You are talking of capacity of two different types of channel.

In one case, the channel inputs and outputs are discrete in time.
At the $i$-th time instant, the received signal is $X_i+N_i$ where
$X_i$ is the received symbol of average energy $E$ and $N_i$ is the
noise (typically modeled as a sequence of i.i.d. $\mathcal N(0,\sigma^2)$
random variables).  The channel capacity of this discrete-time
Gaussian channel $~$is
$$C = \frac{1}{2}\log_2\left ( 1 + \frac{E}{\sigma^2}\right) ~\text{bits per channel use}$$
and so
bits per channel use is the natural metric.  If we are told how far
apart in time the discrete time instants are, e.g. one channel use per
microsecond, then a capacity $C$ bits per channel use can be stated
as bits per second, e.g. $C$ Mbps for our one microsecond
example.

In the second case, the inputs and outputs are continuous-time signals which occupy bandwidth and so the natural measure is bits per second per Hertz.  There
are more complications involved in making the transition from the continuous-time
channel to the discrete model, and in connecting the bandwidth $W$, the received
signal $P$ and the noise spectral density $N_0$ to $E$ and $\sigma^2$
(see here
for some details), but when all
this is done we get Shannon's celebrated formula
$$C = W\cdot \log_2\left ( 1 + \frac{P}{N_0W}\right) ~ \text{bits per second}$$
for the capacity of the additive white Gaussian noise (AWGN) channel of
bandwidth $W$.  This capacity can also be expressed as $C/W$ bits per
second per Hertz.


oversampled coefficient for existing exponential smoothing Say I have an exponential smoothing for certain $\Delta t$, $t_{i+1} = t_i + \Delta t$. In this sampling, I choose a particular $\alpha$ to filter signal $z_i$ like
$$
v_1 = z_0 \\
v_{i+1} = \alpha\:z_i + (1-\alpha)\:v_i
$$

and we have some sort of agreement as to what $\alpha$ to use.

Now, someone comes along and tells me that I am able to get, say, five times more samples in the same $\Delta t$ period. I need to modify the smoothing factor $\alpha$. It seems to me, intuitively, that the new factor, say $\alpha'$, that has the same discount effect for previous values but with five times more iterations should be
$$
(1-\alpha')^5 = 1-\alpha
$$

Is this correct? Is there any principled approaches to justify this?

Note: The formula above seems sensible in the sense that the original asymptotic variance (for a constant model)
$$
\mathbb{V}ar[v] \to \frac{\alpha}{2-\alpha} \sigma^2
$$
becomes
$$
\mathbb{V}ar[v'] \to \frac{\alpha'}{2-\alpha'} \sigma^2
$$
which, under the equation above, is actually smaller (I checked this). This makes sense as more data would give more precision for $v_t$ as estimator of a constant $z_t = k + \epsilon$, where $\epsilon \sim \mathcal{N}(0,\sigma)$.

I think it's instructive to take a look at the impulse responses, taking absolute time into account. I'm assuming the following recursion:

$$y_n=\alpha x_n+(1-\alpha)y_{n-1}\tag{1}$$

The impulse response corresponding to the recursive system described by (1) is

$$h_n=\alpha(1-\alpha)^n$$

Assume now that we have a discrete-time system with sample interval $T$, and another one with sample interval $T^{\prime}=T/L$. Let $\beta=1-\alpha$, and let $\alpha^{\prime}$ and $\beta^{\prime}$ be the constants of the system with sample interval $T^{\prime}$. In order for the two systems to have the same time constant and the same scaling we require (considering a point in time $t_n=nT=nLT^{\prime}$)

$$\alpha\beta^{n}=\alpha^{\prime}{\beta^{\prime}}^{nL}\tag{2}$$

From (2) we see that

$$\alpha=\alpha^{\prime},\quad \beta=\beta^{\prime L}\quad
\text{(i.e. $\beta^{\prime}\neq 1-\alpha^{\prime}$)}$$

In this case the impulse response of the system with the higher sampling rate is simply an interpolated version of the other system's impulse response. Note that the condition on $\beta^{\prime}$ is the one that you came up with intuitively, and it is also the more important of the two because it affects the time constant of the system. The choice $\alpha=\alpha^{\prime}$ only affects the scaling and in this sense is secondary. Choosing $\alpha^{\prime}=1-\beta^{\prime}$ is of course also possible and just changes the scaling.


identify equilibrium region after large transient I am looking at data from a mechanical system under an unsteady load. I'm trying to find the simplest way to identify the portion of the signal once the system reaches its new equilibrium. Here's a simple example.

Any suggestions for a simple and robust approach would be most welcome!

It is difficult to think about robustness without knowing anything about the noise characteristics.

For example, it might be enough to obtain a number of these waveforms (20-50), obtain a box plot of the amplitude at each time instance and use that information to derive a very simple threshold below which the signal enters a region that can be safely characterised as equilibrium.

If there is enough noise to cause the signal to "jump" above the threshold, then you might want to try introducing some form of hysteresis in the Thresholding process. In other words, you might consider "equilibrium" the point when the waveform has dropped below a certain threshold and has stayed there for at least a $u$ amount of time.

If the noise characteristics are such that it is impossible to "tune" such a threshold then a simple 3-5 tap one-dimensional median filter or a very simple low pass filter might be enough, to smooth out the waveform sufficiently and enable you to apply a threshold more reliably.

If the process is such that the waveform is "floating" and a fixed threshold is impossible, then you would be looking at the point when the rate of decay is sufficiently flat. Since you already know that you have an exponentially decaying waveform, you should use this to your advantage. In this case, you can still apply some sort of smoothing pre-processing (e.g. filtering, if required) and then fit the data to an exponentially decaying model. You can then obtain the first order finite difference and set a threshold on that variable for when it drops sufficiently close to zero denoting that the curve is flattening out.

But in your case, you seem to have some sort of undershoot as well (some sort of inertia effect (?)). In extreme cases, the exponential decay model might not be a good assumption in terms of fitting the exact waveform. In that case, it might be better to do a DFT on the settling segment and try to quantify the settling time by known results from automatic control system theory.

Hope this helps.


Why is TV static noise always black and white? Most modern day Cathode Ray Tube (CRT) televisions manufactured after the 1960s (after the introduction of NTSC and PAL standards) supported the circuit-based decoding of colored signals. It is well known that the new color standards was created to permit the new TV sets to be backwards compatible with old black and white broadcasts of the day (among also being religiously backwards compatible with numerous other legacy features). The new color standards added the color information on a higher carrier frequency (but at the same duration of the luminosity). The color information is synchronized after the beginning of each horizontal line and is  known as the colorburst.

It would seem that when you feed noise into a television, the TV should create not only black and white noise but also color noise as there would be color information at each new horizontal line where each frame should be. But this is the not the case as all color TVs still make black and white noise!

Why is this the case?

Here is an example signal of a single horizontal scan.

And here is the resulting picture if all horizontal scans are the same (you get bars!).

The color burst is also an indicator that there is a color signal.

This is for compatibility with black and white signals. No color burst means B&W signal, so only decode the luminance signal (no croma).

No signal, no color burst, so the decoder falls back to B&W mode.

Same idea goes to FM stereo/mono. If there is no 19 kHz subcarrier present, then the FM demodulator falls back to mono.


Gender classification using Periocular Region I have a dataset of the periocular region, I have images of the male and female periocular region but they are only labeled as left and right periocular region, Folder is a mixture of male and female periocular region. One can not judge which one is of male and which one is of female periocular region.

According to biometric classification male and female have different periocular regions, so in short, I have only 2 different types of the periocular region (male and female) in my folder as the dataset.

How I can classify them that system automatically separates the feature of the male and female periocular region and train them as 2 different classes.

You can't if you have nothing to start labelling them with.

As you describe it, there's simply no info available (not even about statistics) about the gender of your reference images, so there's nothing you can train.

What you can do is use a different technique (as you say "biometrics say..." I assume there's a non-ML method of classification) and run it over the dataset, and use the resulting labels to train your neural network. Discussion of the sensibility of that is left to you as an exercise for your report!


How to compute autocorrelation of signal defined by difference equations? I have no experience with difference equations and I want to learn how to compute the following, but I found no resource online. Any help would be greatly appreciated.

Find:

$$\mathbb{E}\left[d[n]d[n + k]\right]$$

Where:

$d[n]$ is a discrete time signal of the form:

$$d[n] = \alpha \,d[n−1]+v[n] \quad (\alpha \in [0;1])$$

$v[n]$ is a sequence of uncorrelated zero-mean unit-variance random variables.

There are different approaches to compute the requested auto-correlations but I would like to provide the simplest that suits to your case provided that.

The difference equation is an LCCDE type indicating an LTI system.
The input to the system $v[n]$ is WSS (wide sense stationary) random process.

Both of which I assume to be observed in your question. Then I would use the following well known relation between the input and output auto-correlations of an (real) LTI system driven by (real) WSS input $v[n]$:

$$ r_{dd}[k] = h[k] \star h[-k] \star r_{vv}[k] $$

where the auto-correlation of WSS $v[n]$ is $r_{vv}[k] = \sigma_v^2 \delta[k] = \delta[k]$ for an uncorrelated, zero mean, unit variance process. And $h[n]$ is the impulse-response of the LTI system signified by the LCCDE.

Then the auto-correlation sequence of the WSS output $d[n]$ is
$$ r_{dd}[k] = h[k] \star h[-k] = \sum_{m=-\infty}^{\infty} h[m]h[-(k-m)]~~~,~~~\text{ for } k = 0,\pm 1, \pm 2,... $$

So assuming that your LCCDE :
$$d[n] - \alpha  ~d[n-1] = v[n] ~~~, ~~~\text{ with } |\alpha| < 1 $$

signifies a causal LTI system, then its impulse response is:

$$ h[n] = (-\alpha)^n u[n] $$

and yields the auto-correlations to be:
$$ r_{dd}[k] = \sum_{m=-\infty}^{\infty} (-\alpha)^m u[m] (-\alpha)^{m-k} u[m-k]~~~,~~~\text{ for } k = 0,\pm 1, \pm 2,... $$

I hope you can proceed the rest.


Extrapolate a 2D array using Fourier Transform I need to extrapolate a given 2D array to a larger domain, keeping the spatial frequency. This is the original field:

(the data file in numpy npz format and a Jupyter notebook to plot it can be found here)

The horizontal size here is 14.864408108 (critical wavelength), and the vertical size is 14.864408108/sqrt(3)

I would like to extrapolate it to a square domain of size, let's say, ten times bigger than the original horizontal size, but fixing the number of grid points to a given number (256x256 for example).

Currently I use numpy.tile to do that:

import numpy as np
from scipy import interpolate

def enlarge(field,mult_x=5,mult_y=8):
    y = np.linspace(2*mult_y*Y.min(),2*mult_y*Y.max(),int(np.round(2*mult_y*field.shape[0]/mult_x)))
    field2 = np.concatenate((field, field[::-1, :]))
    field4 = np.concatenate((field2, field2[:,::-1]), axis=1)
    large = np.tile(field4,(mult_y, mult_x)) # 17 just because it gives a number close to 256*10
    less = large[::mult_x, ::mult_x] # skip every 10 rows/columns
    #print(less.shape)
    f = interpolate.interp1d(y, less,axis=0)
    y_new = np.linspace(2*mult_y*Y.min(),2*mult_y*Y.max(),less.shape[1])
    field = f(y_new)
    return field

And I get this figure:

Can I use Fourier transform to do that?

You can use the Fourier transform for a lot of things – it can be part of a particular step in a scaling algorithm. Still, what you're looking for are scaling algorithms.

In this particular case, your original data size is (74 × 128) – what that means in wavelengths or any physical unit doesn't matter to algorithms – and you want to scale it to (512 × 512).

So what you'd need to get 128 -> 512 would be simply inserting 3 zero rows after every image row and then low-pass filtering the result. A 128-FFT, followed by a zero-padding the 512-128 higher bins and then doing a 512-IFFT would do that – but not like you'd want it to, but with circular effects at the border; what you'd be doing mathematically would be a cyclic convolution with a 1/4 band-wide filter.

Same goes for the 74 -> 512 dimension – but here you'd need to interpolate higher first (to 37·512, to be exact, i.e. by a factor of 37·256) just to decimate down back to 512 again. And then you'd still be stuck with the circular convolution.

What you in fact want to do is apply resampling by a rational rate – the core filter for that could be applied by fast convolution, i.e. with the help of overlapping multiple FFTs, but it's not really the point of the algorithm (just like "using a screwdriver" isn't the central part in "building an airplane").

So, look for any of the existing image scaling algorithms.


Non-linearly additive spectral subtraction? The following pictures describe my situation. Any thoughts on an appropriate algorithm I could use to solve this problem?

It's difficult to suggest an algorithm without knowing more about the nature of these signals, but if it really is as simple as the above situation, then a spectral subtraction approach which uses the reference 'Single Substance' signal to activate a binary mask might do the job for you.


Are there analog re-configurable filters? I know that we can build a chip (or circuit) or use FPGA to have re-configurable digital filters. But how about analog filters, can we have re-configurable analog filters?
thanks

can we have re-configurable analog filters?

Yes. The knob you turn on your grandma's kitchen or living room radio changes the tuning of an oscillator by changing the capacitance of a component. Any stable oscillator is essentially a filter for its resonant frequency.
Also, plenty of other examples: Tunable RC filters with adjustable R; mechanically tunable metal spring filters; tunable microwave cavity resonators with tuning screws, PIN diode-based LC filters... the list goes on.
Basically, take any analog filter topology and make one of the components adjustable, you get a reconfigurable filter.
Is that useful? Well, sometimes it is, if you need to keep your passband adjustable (though that is more commonly done via octave filter banks) or need to adjust your ADC's anti-aliasing filter (though that is more commonly done using high-rate ADCs and digital filtering / sampling methodology).
Generally speaking, you avoid doing complicated things with filters in analog domain as far as possible, really. Analog filters are always a trade-off, and digital filtering hence has replaced them wherever feasible.
Can you get as much freedom as with a digital filter? Practically never. Something like a RRC is pretty trivial in digital, but will need a lot of work to even remotely approximate in analog domain, it will be dependent on temperature, and costly. The closest you'd get to an analog RRC would probably be a bucket-brigade filter, and that's almost a digital filter: discrete in time, but continuous in amplitude.


Zero-phase-shift filters must be symmetric about the origin In reading Rafael Gonzalez Digital Image Processing, second edition, section 4.10.2,
the author writes:

  Notch filters are the most useful of the selective filters. A notch filter rejects (or passes) frequencies in a predefined neighborhood about the center of the frequency rectangle. Zero-phase-shift filters must be symmetric about the origin, so a notch with center at $(u_o, v_o)$ must have a notch at location $-(u_o, v_o)$.

Two questions:

1) Why does a zero phase shift filter have to be symmetric about the origin (with real coefficients I am assuming)?

2) What does the author mean here when he says a "notch"

Any insights appreciated.

Consider the Inverse Fourier Transform as the impulse response of the filter. In the case of an ideal brick-wall filter with zero-phase, the impulse response would be a Sinc function centered about zero in time (given by the Inverse Fourier Transform).

This is non-causal (time domain response is greater than 0 for $t<0$), but in order to be zero-phase as shown, the time domain response must be complex conjugate symmetric.  It need not be real; in order for the frequency response to be real, the time response real components must be even and the imaginary components must be odd (complex conjugate symmetric).

But to visualize the above plot further and what "zero-phase" means, a delay in time is a linear phase in frequency, thus if we shift the impulse response to the right (as done prior to truncation for a causal solution), the magnitude response is unchanged but the phase in no longer "zero-phase". If we continued to see the next plot I have in this series, which shows the impulse response truncated by the dashed window in this plot below, we would then see that the magnitude response of the filter would also be distorted, but that is beyond the question being asked here.

I summarize the properties of odd and even functions as well as other universal Fourier Transform properties in the figure below. The time and frequency domains are interchangeable for these properties; for example, if a function is periodic in one domain, it is discrete in the other. If a function has only real components in one domain (zero-phase), then all real components in the other domain must be an even function about 0, and all imaginary components must be odd about zero (inverse sign). This results in complex conjugate symmetry.

The relationship of the Hilbert Transform mentioned for the causal and anticausal functions applies to minimum phase systems.

What is interesting is how these tie together and why I choose to describe the above relationships in terms of even and odd functions rather than simply saying the other domain will be complex conjugate symmetric. This way we can see how a causal function in time (or a single sided function in frequency) is the sum of even and odd functions, and therefore MUST be complex in the other domain given that we need both real and imaginary components to create such even and odd functions.

Regarding the OP’s second question: a notch filter usually implies a very narrow band of rejection to remove certain frequency components. For more details see this post here: Transfer function of second order notch filter


Meaning of wavelet and scaling coefficients What is the meaning of wavelet coefficients and scaling coefficients? E.g. for a sequence I obtained the following wavelet coeffients.

How am I supposed to interpret them? I used wavelets package in R. The sequence concers price time series of 21 days. I used the Haar wavelet.

Hope this helps:...

The fast wavelet transform (FWT) is a mathematical transform which results in a series of orthogonal coefficients of varying resolution.

Assume an input signal or time series $S$ of length $n=2^{j}$.  The length of $S$ must be equal to some power of 2, that is, $n$ equal to 8, 16, 32, 64, 128,  256, 512, 1024, etc.  To transform a signal into wavelets, at each resolution the signal is convolved using a series of scaling coefficients, $a_0$, $a_1$, $a_2$, $\ldots$ representing a low pass filter and a sequence of $d_0$, $d_1$, $d_2$, $\ldots$ of a high-pass filter.

Scaling function coefficients, $\phi$.      At the $j-1$th resolution, the convolution is taken over the signal's original length, $2^j$, in the form
\begin{equation}
\phi_k =  a_{-1}s_{2k+1} + a_0 s_{2k},
\end{equation}
where $k=0,1,2,\ldots,2^{j-1}-1$ are the indices of each scaling function, $s_0, s_1, \ldots, s_{2^{j}-1}$ are the original input signal elements, and $a_{-1}$ and $a_0$ are constants from the Haar scaling function.  An example of expanding the notation for the $j-1$ resolution is
\begin{equation}
\begin{split}
\phi_0 &= a_{-1}s_1 + a_0 s_0\\
\phi_1 &= a_{-1}s_3 + a_0 s_2\\
\vdots\\
\phi_{2^{j-1}-1} &= a_{-1}s_{(2^j-1)} + a_0 s_{(2^j-2)}.\\
\end{split}
\end{equation}
Consider a signal $S$ of length $n=1024=2^{10}$ $(j=10)$, and assume Haar scaling function values of $a_{-1}=a_0=0.5$.  Although the length of $S$ is 1024, the notation used to represent the series of values is $s_0$, $s_1$, $s_2$, $\ldots$, $s_{1023}$.   The first convolution of $S$ will result in 512 ($2^{10-1}=2^9$) scaling function coefficients, which are obtained using the relationship
\begin{equation}
\phi_k = 0.5 s_{2k+1} + 0.5 s_{2k}
\end{equation}
where $k=1,2,\ldots,511$.  This translates to
\begin{equation}
\begin{split}
\phi_0 &= 0.5 s_{1} + 0.5 s_{0}\\
\phi_1 &= 0.5 s_{3} + 0.5 s_{2}\\
\vdots\\
\phi_{511} &= 0.5 s_{1023} + 0.5 s_{1022}\\
\end{split}
\end{equation}
It is easily noticed that the scaling function coefficients are the average of each neighboring pair of signal elements.

Wavelet coefficients, $\psi$.  The wavelet coefficient is essentially based on the difference between each neighboring pair of signal elements.  At the $j-1$ resolution, these are
\begin{equation}
\begin{split}
\psi_0 &= -d_{-1}s_1 + d_0 s_0\\
\psi_1 &= -d_{-1}s_3 + d_0 s_2\\
\vdots\\
\psi_{2^{j-1}-1} &= -d_{-1}s_{(2^{j-1}-1)} + d_0 s_{(2^{j-1}-2)}\\
\end{split}
\end{equation}
where $d_{-1}$ and $d_0$ are both 0.5 based on the Haar wavelet.

For our first convolution of signal $S$ of length 1024 at resolution $j-1$, the 512 wavelet coefficients are
\begin{equation}
\begin{split}
\psi_0 &= -0.5 s_{1} + 0.5 s_{0}\\
\psi_1 &= -0.5 s_{3} + 0.5 s_{2}\\
\vdots\\
\psi_{511} &= -0.5 s_{1023} + 0.5 s_{1022}\\
\end{split}
\end{equation}


Phase spectrum: $\mathrm{arg}\mathcal{F}(\omega)$ or $\mathrm{-arg}\mathcal{F}(\omega)$? Do we need minus or not?
I need to extract phase spectrum from that thing

$$f(x) = \frac{1}{2\pi}\int_{-\infty}^{\infty} e^{i\omega x}d\omega\int_{-\infty}^{\infty}f(t)e^{-i\omega t}dt$$

I can do it by taking Fourier Transform of $f(x)$, right? And phase spectrum will be:

$$S(\xi)=\int_{-\infty}^{\infty}f(x)e^{-i\xi x}dx $$

$$\Phi (\xi)=-\mathrm{arg}\,S(\xi)$$

I found that in book, but don't understand do we need minus
 before arg or not? I'm looking on different plots and some sources plotting with minus and some without. Help me please

Both definitions are common. However, if you talk about the phase, then what is usually meant is the (frequency dependent) phase angle $\phi(\omega)$ of the Fourier transform $F(\omega)$:

$$F(\omega)=|F(\omega)|e^{j\phi(\omega)}\tag{1}$$

The negative phase $\beta(\omega)=-\phi(\omega)$ is often called the phase shift or phase lag. It is commonly used together with the attenuation $\alpha(\omega)$, which is the negative natural logarithm of $|F(\omega)|$:

$$\alpha(\omega)=-\log (|F(\omega)|)\tag{2}$$

$$F(\omega)=e^{-(\alpha(\omega)+j\beta(\omega))}\tag{3}$$


What happens if an harmonic in a band-limited signal at Nyquist frequency is added to a 90 degrees out-of-phase replica of it? I'm interested about aliasing in digital audio and I wonder if aliasing can be produced by simple mixing of band-limited signals.
As I know a band limited signal can contain frequencies up to sample rate/2 frequency, the so-called Nyquist frequency. Is this right?
If that is so, what will happen if an harmonic in a band-limited signal at Nyquist frequency is added to a 90 degrees out-of-phase replica of it?
Will this create an harmonic with double frequency and could such an addition create aliasing?

No.  Adding two waves at the same frequency results in another wave of the same frequency.  Depending on the phases, there can be anywhere from complete constructive interference to maximum destructive interference.  If the amplitudes are the same there can be complete cancellation.  This applies to the Nyquist frequency too.

See the math from my answer to FFT Frequency bin relationship.

Hope this helps.

Ced


can i decompose image deblurring model? i am working in digital image restoration field, recently i have studied image degradation model which is like that.
g(x,y) = h(x,y)*f(x,y) + n(x,y)

where

g(x,y) is the degraded image

f(x,y) is the image that we wish to capture or its is latent image

h(x,y) is image degrading function, and * is convolution operator

n(x,y) is noise

my question is that, can i change this particular model like that.

g(x,y) = h(x,y)*{f(x,y)+n(x,y)}

In principle you can do that. You should, however, note that the two noise signals (both of which you called $n(x,y)$) are of course not identical. The two models are

$$g(x,y)=h(x,y)*f(x,y)+n(x,y)\\
g(x,y)=h(x,y)*\left[ f(x,y)+z(x,y)\right]$$

where the two noise components $n(x,y)$ and $z(x,y)$ have different properties. So you have to be careful with methods that assume certain noise properties. E.g. if it is assumed that $n(x,y)$ is white, then $z(x,y)$ is not white anymore and this has to be taken into account in the noise removal algorithm. The basic question you have to ask yourself is what you gain by changing the model.


Matching features points of rotate image using SURF I'm trying to match SURF points between 2 images using SURF with opencv. One is the rotation of the other.
The problem is that the few matching it can found are wrong.

image of matching points

I mainly took my code from opencv http://docs.opencv.org/doc/tutorials/features2d/feature_flann_matcher/feature_flann_matcher.html#feature-flann-matcher

Here is the code:

string imageName1="test_right_rotate.jpg";
string imageName2="test_right.jpg";
Mat image1 = imread( imageName1, 1 );
Mat image2 = imread( imageName2, 1 );
int minHessian = 500;
SurfFeatureDetector detector( minHessian );
std::vector<KeyPoint> keypoints_1, keypoints_2;

detector->detect( image1, keypoints_1 );
detector->detect( image2, keypoints_2 );

SurfDescriptorExtractor extractor;
Mat descriptors_1, descriptors_2;

extractor.compute( image1, keypoints_1, descriptors_1 );
extractor.compute( image2, keypoints_2, descriptors_2 );

FlannBasedMatcher matcher;
std::vector<vector<DMatch > > matches;
matcher.knnMatch(descriptors_1,descriptors_2, matches, 2);
std::vector< DMatch > good_matches;
for(int i = 0; i < min(descriptors_1.rows-1,(int) matches.size()); i++) //THIS LOOP IS SENSITIVE TO SEGFAULTS
{
    if((matches[i][0].distance < 0.6*(matches[i][1].distance)) && ((int) matches[i].size()<=2 && (int) matches[i].size()>0))
    {
        good_matches.push_back(matches[i][0]);
    }
}
Mat img_matches;
drawMatches( image1, keypoints_1, image2, keypoints_2, good_matches, img_matches, Scalar::all(-1), Scalar::all(-1), vector<char>(), DrawMatchesFlags::NOT_DRAW_SINGLE_POINTS );

It works well for scaled images and translated images but fails for rotated images.

EDIT:

I don't know what is wrong but instead of using SURF, I used ORB as detector and FREAK as extractor and it works a lot better.
You can find the code on this SO : https://stackoverflow.com/questions/12491022/opencv-freak-fast-retina-keypoint-descriptor

SURF Features are not completely rotation invariant.

In the paper, they try to make SURF more robust to rotations by following the following methodology.

In order to improve the performance for rotation, the descriptors are based on multi-orientations. The many-to-many tentative correspondences are determined with a maximum distance. Hough transform is used to reject the mismatches and the affine parameters are computed with a least-squares solution.

The image matching algorithm shows a better performance for image rotation than the standard SURF and it succeeds in matching the image including repetitive patterns which will deteriorate the distinctiveness of feature descriptors.


Finding the Transfer Function of a Multiplicative Distortion in MATLAB I have the frequency response of two signals shown above. One is uncontaminated, and the other is the same signal contaminated by a multiplicative distortion with two poles and two zeros. I need to find/estimate the transfer function values of the contaminating filter in order to equalize the contamination. How do I accomplish this in MATLAB? Thanks.

Basically the transfer function is given by (Continuous time case):
$$ H \left( j \omega \right) = K \frac{ \left( j \omega - {q}_{1} \right) \left( j \omega - {q}_{2} \right) }{ \left( j \omega - {p}_{1} \right) \left( j \omega - {p}_{2} \right) } $$
Now, given your frequency data sampled an $ N $ points $ {\left\{ {\omega}_{i} \right\}}_{i = 1}^{N} $ we need to use the discrete-time form of the transfer function over the Z domain (Z Transform).
Hence the model becomes:
$$ H \left( z \right) {\mid}_{z = {e}^{j \omega}} = K \frac{ \left( z - {q}_{1} \right) \left( z - {q}_{2} \right) }{ \left( z - {p}_{1} \right) \left( z - {p}_{2} \right) }{\Bigg|}_{z = {e}^{j \omega}} = K \frac{ \left( {e}^{j \omega} - {q}_{1} \right) \left( {e}^{j \omega} - {q}_{2} \right) }{ \left( {e}^{j \omega} - {p}_{1} \right) \left( {e}^{j \omega} - {p}_{2} \right) } $$
you may set a system of $ N $ complex equations and use non linear solver to find the 5 parameters ($ K $, $ {q}_{1} $, $ {q}_{2} $, $ {p}_{1} $, $ {p}_{2} $)  you need.
In case the measurements are noisy you may use Least Squares.
Remark
The above solves the problem from the data of @jmcd in frequency domain. In practice, using the data domain one could use simpler System Identification (See System Identification: An Overview) methods to get the equivalent ARMA model.
Regarding Frequency Domain based methods for System Identification have a look at Rik Pintelon - System Identification: A Frequency Domain Approach.


PSK vs FSK for noise immunity? Which method provides better noise immunity?

I know you can transmit more data via psk, but what about which method is less susceptible to noise?

As you say in your comment, PSK signals are more bandwidth efficient, but FSK signals have good noise rejection for their data rates.  It is thus difficult to compare them on even terms.
The traditional way of measuring noise immunity is measuring/calculating the probability of a bit error given the bit energy vs. noise energy (Eb/No).  This is a good way of making everything fair
The following plot uses this paper for the FSK analysis, and this page for the BPSK analysis.

As you can see, BPSK wins quite handily.  I won't attempt to prove it, but QPSK has the exact same BER curve as BPSK.  You can intuitively see that if you think about how the real and imaginary portions of the symbol are orthogonal to each other, so you are essentially just sending two orthogonal BPSK signals at the same time.  The BER curves of the higher PSK types (e.g. 8-PSK), though, do get worse because the symbols are closer together, so more energy is needed to distinguish them in the presence of noise.
So why does anyone use FSK?  Mostly because it's pretty simple to implement, particularly the non-coherent receivers.


Kalman filter, defining the measurement model I would like to implement a Kalman filter to estimate the velocity and position of an object. I have an accelerometer, therefore the acceleration is known. The approach is same as:

Kalman filter for position and velocity: introducing speed estimates

however, I don't understand well, how I should define the measurement model:

$\hat{X}(k+1) = \hat{X}(k) + K * ( y(k) - H*\hat{X}(k) );$

K is the kalman gain and is calculated in every step, but What is my "y" here, i know it is called the observations, but how i have my observations ? should i calculate the:

$v(k+1) = a(k)*t + v(k);$

$x(k+1) = 0.5*a(k)*t^2 + v(k)*t +x(k);$

Thank you!

Your observation are the $y_k$ your observation model is defined by $Y_k=HX_k+\eta_k$, where $\eta_k$ if a sample from a normal distribution.

If you want your kalman filter to estimate position, velocity and acceleration your state vector is : $X_k=[x_k,\dot{x}_k,\ddot{x}_k]^T$ and therefore your observation matrix is : $H=[0,0,1]$.

As stated by Rhei, your observations are what you get from the system at each time step, i.e. the measured acceleration.


Sum of two random noises Let $u_1$ be a random noise in $\mathbb R$. The PSD of $u_1$ is constant in the interval $[f_{1,1}; f_{1,2}]$ and is noted $p_1$.
Let $u_2$ be a random noise in $\mathbb R$. The PSD of $u_2$ is constant in the interval $[f_{2,1}; f_{2,2}]$ and is noted $p_2$.
The intervals $[f_{1,1}; f_{1,2}]$ and $[f_{2,1}; f_{2,2}]$ are disjoint.
Is the PSD of $u_1+u_2$ computable ? Does it only depends on $p_1$ and $p_2$?

Is the PSD of u1 and u2 computable ? Does it only depends on p1 and p2 ?

Yes. The Fourier Transform is a linear operation so the spectrum of the sum is simply the sum of the spectra. The PSD is just magnitude squared of the Fourrier Transform but since the frequency intervals are disjoint the squaring doesn't create any type of interaction between the two original spectral. So it's simply
$$P_{sum}(f) = \begin{cases}
 & p_1\text { if } f_{11} < f < f_{12} \\
 & p_2 \text{ if }  f_{21} < f < f_{22} \\
 & 0 \text{ otherwise }
\end{cases} $$


Polyphase filter notation Hi I am a little confused on what the notation of the following statement means.

$$ H_{k}(z)= H(W_{4}^{k} z), k = 0,...,3$$

It comes from a question in which I have designed a FIR low-pass filter $H(z)$ and my goal is to implement a DFT filter bank scheme like this:

Exchange $P(z)$ for $H(z)$ and k would correspond to the subscript of P and M in this case is equal to 4

I guess I am confused on how to find $H_{k}(z)$ or what exactly a polyphase filter is

Hk are modulations of the low pass filter (band pass instead of low pass).

$$
(W_{4}^{k}) = e^{-2j\pi k /4}.
$$
For $z = e^{j\omega}$: $$H(W_{4}^{k} z) = H(e^{j(\omega-2\pi k/4)})$$
This means that the filters $H_k$ are shifted in frequency- these are the band pass filters you want to get using your filter bank.
For $k=0$ $H_0$ will pass the frequencies $[\frac{-\pi}{4} \frac{\pi}{4}]$, for $k=1,[\frac{\pi}{4} \frac{3\pi}{4}]$ etc..

In the DFT filter bank scheme,
$y_0[m]$ are the outputs from $H_0$ , $y_1[n]$ are the outputs from $H_1$ and so on.


How to find f given that f(some_image) = processed_image I have a before/after image pair with pixel transformations only (for example, grayscale, or a hue shift. No blur, border, sharpen, etc.). I'm trying to find a function such that for every pixel (r, g, b) in the original image, f(r, g, b) = (r', g', b'), (r', g', b') being the pixel in the "after" image.

I was thinking about a polyfit function, but I couldn't find one that works in 3D.

I'd appreciate any help with this, or another method of achieving what I need.

Fit 3 polynoms (or any other model) using the individual channels of the output image: r' as a function of (r, g, b), g' as a function of (r, g, b), b' as a function of (r, g, b). I would also suggest not including in the data set all pixels with a 0 or 255 value (if your image is 8-bit) - these values might have been clipped.


Resonant peak frequencies for phaser effect I have a conventional phaser setup with allpass filters and feedback.

I can calculate the peaks frequencies of the phaser when the allpass delays are set to 1 (then they work as simple one sample delays). Adding two allpass delays will give an additional peak, and the peak frequencies for peaks i = 1 to N is found by:

(i / (N + 0.5) ) * SampleRate / 2

where N is the number of peaks (and the number of allpass-delays divided by 2).
So if you have a 6-stage phaser at 44100, the 3 peaks will be at:

(1 / (3 + 0.5) ) * 44100 / 2 = 6300
(2 / (3 + 0.5) ) * 44100 / 2 = 12600
(3 / (3 + 0.5) ) * 44100 / 2 = 18900

However, these frequencies start warping when the allpass-delays are not set to 1:

How can I calculate the peak frequencies when the allpass-delays are not set to 1?

I refer to the following structure of a phaser (from this article):

Let's call the depth factor $a$, and the feedback factor $f$ ($|f|<1$).
I assume that there are $2N$ first order allpass filters, each with transfer function

$$A_0(z)=\frac{z^{-1}-b}{1-bz^{-1}},\quad |b|<1\tag{1}$$

The allpass frequency response is obtained by evaluating (1) on the unit circle $z=e^{j\omega}$:

$$A_0(e^{j\omega})=\frac{e^{-j\omega}-b}{1-be^{-j\omega}}=
\frac{e^{-j\omega}-2b+b^2e^{j\omega}}{|1-be^{-j\omega}|^2}=\\
=\frac{(1+b^2)\cos(\omega)-2b+j(b^2-1)\sin(\omega)}{|1-be^{-j\omega}|^2}\tag{2}$$

And since $A_0(z)$ is an allpass function

$$A_0(e^{j\omega})=e^{j\phi(\omega)}\tag{3}$$

From (2) and (3) the phase $\phi(\omega)$ of the allpass satisfies

$$\tan(\phi(\omega))=\frac{(b^2-1)\sin(\omega)}{(1+b^2)\cos(\omega)-2b}\tag{4}$$

If you cascade $2N$ of these allpass filters, the frequency response of this cascade is

$$[A_0(e^{j\omega})]^{2N}=e^{j2N\phi(\omega)}$$

Let us first consider a system without feedback, i.e. $f=0$. The total frequency response of such a system is

$$H(e^{j\omega})=1+ae^{j2N\phi(\omega)}=
1+a\cos(2N\phi(\omega))+ja\sin(2N\phi(\omega))$$

and its squared magnitude is

$$|H(e^{j\omega})|^2=1+2a\cos(2N\phi(\omega))+a^2\tag{5}$$

From (5) it is obvious that the locations of the maxima of $|H(e^{j\omega})|^2$ are independent of the depth paramter $a$, and that they only depend on $N$ and on the allpass phase response $\phi(\omega)$ (i.e. on the paramter $b$). The maxima occur where the argument of the cosine in (5) equals $2\pi k$, $k=0,1,\ldots,N$, which is equivalent to

$$\phi(\omega_k)=\frac{k\pi}{N},\quad k=0,1,\ldots,N$$

where $\omega_k$ are the peak frequencies. With (4) this condition becomes

$$\frac{(b^2-1)\sin(\omega_k)}{(1+b^2)\cos(\omega_k)-2b}=\tan\left(\frac{k\pi}{N}\right),\quad k=0,1,\ldots,N\tag{6}$$

Equation (6) can be used to determine the peak frequencies $\omega_k$. Note that this needs to be done numerically since this equation cannot be solved for $\omega_k$ analytically. There are just some special indices $k$ for which an explicit solution is possible:

$$\begin{array}{rcl}\omega_0&=&0\\
\omega_{N/2}&=&\arccos\left(\frac{2b}{1+b^2}\right)\quad\textrm{if } N \textrm{ is even}\\
\omega_N&=&\pi
\end{array}$$

The following figure is an example for $N=4$, $a=0.7$, and $b=0.6$ which shows $|H(e^{j\omega})|$
and the locations of the maxima (shown by 'x') which I've computed using equation (4).

If there is feedback ($f\neq 0$) then the frequency response looks more complicated:

$$H(e^{j\omega})=\frac{1+(a-f)e^{j2N\phi(\omega)}}{1-fe^{j2N\phi(\omega)}}$$

and its squared magnitude is

$$|H(e^{j\omega})|^2=\frac{1+2(a-f)\cos(2N\phi(\omega))+(a-f)^2}{1-2f\cos(2N\phi(\omega))+f^2}$$

Luckily it turns out that feedback doesn't change the location of the peaks, it just makes them 'peakier'. In the following figure the system the frequency response shown above is compared to the response with feedback ($f=0.5$):


Design of butterworth filter I was trying to design a butterworth filter of arbitrary order. Then i came across this presentation:

Lecture 9: Poles, Zeros & Filters, by Peter Cheung, Imperial College London

I understand most of it, except for 2 things,

Now, we use the fact that $-1=e^{j\pi(2k-1)}$ and $j=e^{\frac{j\pi}{2}}$. I have no idea what's he saying here.
And does he take a normalized cut-off frequency? If yes, how do i change it to arbitrary cut-off frequency.

Now, we use the fact that $-1=e^{j\pi(2k-1)}$ and $j=e^{\frac{j\pi}{2}}$. I have no idea what's he saying here.

Well, he's simply stating mathematical truths:

$f(x)=e^{jx}$ is periodic with $2\pi$, and $e^{j\pi(0-1)}=e^{-j\pi}=-1$. If that's not clear, draw a unit cycle in the complex plane and remember Euler's formula.
$j=e^{\frac{j\pi}{2}}$: same business. $\frac\pi2$ is a quarter cycle, and thus you end up at the top point of the unit circle.

You really need to be familiar with the properties of the complex exponential function, or else you will have a very hard time understanding pole/zero diagrams in general.

  And does he take a normalized cut-off frequency?

Yes, he does (very explicit on slide 9).

  If yes, how do i change it to arbitrary cut-off frequency.

Um, same slide, explicitly:

    We can design filters for any other cut-off frequency by substituting $s$ by $\frac s{\omega_c} $.


Variable Step Size LMS vs Leaky LMS Adaptive Filter Algorithm What is the advantage of Variable step size LMS over Leaky-LMS adaptive filter algorithm? Which one has a better performance?

Variable step size LMS is generally used to improve speed of convergence or decrease steady state error.
Leaky adaptation is used to combat problems like potential instability of the filter in a finite-precision implementation. It is closely related to L2 norm regularization technique and results in continuous down scaling of filter coefficients (hence smoothing the extracted filter).
So comparing these two techniques is not that meaningful. In fact you can easily use both techniques together to achieve the above advantages.


How to remove color cast from HSV information I have hue, saturation and value of a color (HSV color-space). The color has to be set on a device that displays them incorrectly. If I specify white it displays light magenta. However, if you specify some green it displays the best white it can (let's say: gray). How do I remove this color cast? Actually: How do I create a color cast in the HSV numbers to compensate for the wrong display?

My approach was the following. But it lead to a skip of some magenta colors, when I ran through the whole rainbow with a saturation of, say, 0.8 (interval 0 to 1):

Shift magenta to the beginning of the hue scale (60 degrees forward) and thus the hue of the current color
Calculate the distance to move in the direction of green. With a saturation of 1 that means to move completely to green (180 degrees).
This is the formula: distance = (180 / 360 - shifted_hue) * (1 - saturation)
Then add the distance to the shifted hue (possibly resulting in a subtraction).
Shift the hue scale back by 60 degrees.

As said, this is leading to an ugly skip of some magenta colors, because the lower the saturation the bigger the gap at magenta.

Do you know of a better algorithm?

It is not necessary that your device's color mismatch can be corrected with a simple shifting of the Hue component (except of course if you are absolutely sure that this is the case with your equipment).

In general, to correct the response of the device you somehow need to recover the mapping that it already performs between reference and observed colors and then see how it can be corrected by some curve (monotonic or not).

For this purpose you will need a colour chart (and ideally, an instrument such as this one).

(Please note that it is best that this correction is applied at the system level, i.e. your application will still be sending some "color demand" normally, without correcting it but that demand would undergo some correction prior to display automatically by the system.)

At the very basic level, you can send a simple reference image to your display and then try to adjust its R,G,B settings (if there are any) so that the colours seem to display properly (according to your eyes that is)

If you have a system that is not capable of applying such profiles but you still must apply some color-correction you could try matching the relative distances of the colours with a known target (or image in your case) and then (if required) try to correct towards absolute colours by shifting (with trial and error).

Take an evenly lit image of the colour chart with a good digital camera, color correct it (to reduce any effects of the lighting or sensor used to acquire your "reference" image) and feed it to your display. Take an image of the display (ideally in the same lighting conditions), apply the same color correction with the profile from the first step (and for the same reasons) and then examine the relative color correspondence between any two colors from the color chart image and two colors from the color chart display image to recover the mapping. (A piece of software that might be of help in this process is LPROF)


find a line go through the center mass and that line must be inside object area I have an object (orange color object) and I know its centroid/center mass C. I want program to automatically find another point ( can be point 1, point 2, point 3 or point 4 and so on) such that the line connect point C and one of those points (so as green line in figure 2) must stay inside the object. Moreover, the green line length must be as long as possible, or at least 3/4 of the maximum length. The yellow point in figure bottom is not qualified because the violet line which connect it to the center mass C has a part which is not inside the object.

Any idea? I am doing it on Matlab (image processing)

How about:

Separate the foreground from the background and keep a list of indices or a logical mask from the background
Convert the image to a (rho, theta) coordinate space, meaning that each pixel is encoded by the angle and length of vector from the origin to itself. You can scan along this vector to see if you encounter any background pixels. Have a look at the so-called Trace Transform
Alternatively you can stay in the original coordinate domain and use Bresenham's line algorithm to get all pixels belonging to a scan-line between a point and the center and then comparing them against the list of background labeled pixels.

Perhaps you can optimize a bit by starting on the edge of the foreground shape, working your way inward to the center, discarding any scan-line the moment it encounters a background labeled pixel.


Is this a "valid" CDMA modulation? I am banging my head against (specific aspects of) CDMA for a few days. At the basic level I am confused about the interrelation between implementations that are modeled as mixing and fully digital implementations. Example:

Mixing could use the the matrix $[ 1, 1 ; 1, -1]$ (appears in the example of https://en.wikipedia.org/wiki/Code-division_multiple_access#Example and https://en.wikipedia.org/wiki/Walsh_matrix). These matrices have elements -1 and 1 and flip signs and can be interpreted as mixing/multiplication. Makes sense
Other sources show these matrices as $[ 0, 0 ; 0, 1]$ (e.g. https://www.tutorialspoint.com/cdma/cdma_techniques.htm). My understanding is that these setups assume bit streams (as opposed to more general sample streams) and for each data "1" they would just send the normal PN sequence (e.g. $[ 0, 1 ]$) and for each data "0" they would send the 1's complement $[ 1, 0]$. Makes sense too.

Now I found a paper with yet another setup:

In the example station one has code $01$ and station two has code $10$. Station one sends $111001$ and station two sends $100111$. After CDMA code injection, station one transmits $010101000001$ and station two transmits $100000101010$. It should be noted that if these two transmissions are summed in the shared channel, the result is still a two-level waveform (=binary): $110101101011$! Demodulation is performed again by multiplying (=ANDing) the chip code and then OR-ing each bit. For example, station two would multiply/AND the received signal $110101101011$, and it with $10$ to get $100000101010$. OR-ing two adjacent bits gives the original sequence $100111$.

Questions:

Is #3 a valid CDMA modulation/demodulation scheme?
How does it relate to the setups #1 and #2 ?
How do setup #1 and #2 relate? (I would expect multiplication in $GF(2)$ is AND and addition is XOR ... but 2. is different)
How would setup #3 extend to code lengths $\gt 2$?
Why would someone use setup #3 and call it "CDMA" (instead of setup #2)? Do I miss the relationship to #2?

One one hand #3 would provide a more logical match with #1 because AND is a multiplication in $GF(2)$. However, addition would correspond to XOR instead of OR and I don't know where OR comes from.

On the other hand, I cannot make setup #3 work with codes that resemble normal Walsh matrices (e.g. $[ 1, 1 ; 0, 1]$. All that works are codes like $[ 1, 0, 0, 0 ; 0, 1, 0, 0 ; 0, 0, 1, 0; 0, 0, 0, 1]$. Yes, they are orthogonal but somewhat miss the point of CDMA: code spreading. They implement more time division.

I think you are trying to mix orthogonality for $GF(2)$ and $R^n$. For the sequence $010101000001$, if station 2 had tried to demodulate it, he will have to take groups of 2 bits, multiple (AND) and add(XOR) them, to decode whether bit 0 or 1. Station 2 would have got all 0s because $10$ and $01$ is 'orthogonal' for station 1 and 2. So yes, #3 is indeed a valid CDMA scheme in $GF(2)$. I stress, it is a valid CDMA scheme only in $GF(2)$. It is not related to #1 and #2 (see below why).

Setups #1 and #2 are valid CDMA schemes in $R^n$ (real number vectors) because orthogonality is validated by $x_1^Tx_2$. For #2, you cannot validate orthogonality by dot product of vectors containing 0 and 1. Because for $W_4$, rows 3 and 4 are not orthogonal. So, $0$ has to be represented by $-1$. Even though logically we switch $0$ and $1$ as NOT operator while forming Walsh matrices, orthogonality is not valid in $GF(2)$.

  On the other hand, I cannot make setup #3 work with codes that resemble normal Walsh matrices

$[1,1;0,1]$ would not correspond to CDMA because whether in $GF(2)$ or $R^2$, they are not orthogonal. $[1,0,0,0;0,1,0,0;0,0,1,0;0,0,0,1]$ is orthogonal in $GF(2)$ as well as in $R^2$ which is same as your multiplexing in time but it is still technically a CDMA scheme (if you transmit each bit at 4x higher rate).


Is a leaky integrator the same thing as a low pass filter? The equation governing a leaky integrator (according to Wikipedia at least) is

$\frac{d\mathcal{O}}{dt} + A\mathcal{O}(t) = \mathcal{I}(t)$.

Is a continuous-time leaky integrator thus the same thing as a low pass filter with time-constant $A$, up to some scaling of the input?

A so-called leaky integrator is a first-order filter with feedback. Let's find its transfer function, assuming that the input is $x(t)$ and the output $y(t)$:

$$
\frac{dy(t)}{dt} + Ay(t) = x(t)
$$

$$
\mathcal{L}\left\{\frac{dy(t)}{dt} + Ay(t)\right\} = \mathcal{L}\left\{x(t)\right\}
$$

where $\mathcal{L}$ denotes application of the Laplace transform. Moving forward:

$$
sY(s) + AY(s) = X(s)
$$

$$
H(s) = \frac{Y(s)}{X(s)} = \frac{1}{s + A}
$$

(taking advantage of the Laplace transform's property that $\frac{dy(t)}{dt} \Leftrightarrow sY(s)$, assuming that $y(0) = 0$).

This system, with transfer function $H(s)$, has a single pole at $s = -A$. Remember that its frequency response at frequency $\omega$ can be found by letting $s=j\omega$:

$$
H(j\omega) = \frac{1}{j\omega + A}
$$

To get a rough view of this response, first let $\omega \to 0$:

$$
\lim_{\omega \to 0} H(\omega) = \frac{1}{A}
$$

So the system's DC gain is inversely proportional to the feedback factor $A$. Next, let $w \to \infty$:

$$
\lim_{\omega \to \infty} H(\omega) = 0
$$

The system's frequency response therefore goes to zero for high frequencies. This follows the rough prototype of a lowpass filter. To answer your other question with respect to its time constant, it's worth checking out the system's time-domain response. Its impulse response can be found by inverse-transforming the transfer function:

$$
H(s) = \frac{1}{s+A} \Leftrightarrow e^{-At}u(t) = h(t)
$$

where $u(t)$ is the Heaviside step function. This is a very common transform that can often be found in tables of Laplace transforms. This impulse response is an exponential decay function, which is usually written in the following format:

$$
h(t) = e^{-\frac{t}{\tau}}u(t)
$$

where $\tau$ is defined to be the function's time constant. So, in your example, the system's time constant is $\tau = \frac{1}{A}$.


How to solve ADMM for TV Minimization Problem For Different Sizes $A$ and $x$ in $Ax=b$ I have matrix $A$ that is $(M \times M)$ square matrix, $x$ $(M \times N)$ matrix, $b$ is $(M \times N)$ matrix. Knowing $A$ and $b$ I would like to get the $x$ from the equation $Ax=b$. $N=p \times q$, so consider $x$ as an $M$ number of $p \times q$ pixel images. These images are sparse when we take total variational $TV(x)$. I would like to solve the following minimization problem for $x$;

Define an operator $g(x)=\parallel x^T\parallel_2,$ or $g(x)=\parallel x^T\parallel_1,$ that maps $M \times N  \mapsto  N \times 1$.

$$\min  \frac{1}{2}\parallel Ax-b\parallel^2_2+k\parallel z\parallel_1,$$ subject to $Fg(x)-z=0$.

$F$ is the difference matrix to take numerical gradient of pixels.

Unlike in the deblurring+denoising problem, my matrix $A$ and $x$ are in different sizes.

The ADMM (Alternating Direction Method of Multipliers) solution to my minimization problem is given as:

$$x^{k+1}=(A^TA+pF^TF)^{-1}(A^Tb+pF^T(z^k-u^k))$$

$$z^{k+1}=S_{t/p}(Fx^{k+1}+u^{k})$$

$$u^{k+1}=u^k+Fx^{k+1}-z^{k+1}$$

When calculating $x^{k+1}$ the matrix dimensions are not compatible in my case. Gradient matris should be applied to total number of pixels $N$, $F$ matrix should have $N$ columns. So the term $(A^TA+pF^TF)$ is not proper way to add 2 matrices.

How can I overcome this and find a solution to my problem?

The Error in the Model

The problem is in the dimensions of the Linear Operator $ A $ in your model compared to the Data Matrix $ X $. The number of columns of the matrix $ A $ must match the number of pixels in each column of $ X $ (Each image). While in your case it matches the number of images.

Matrix Form of 2D Linear Operator

Let's try to build the model correctly. Assume our data (Images) is given by the set of 2D matrices $ {\left\{ {X}_{i} \right\}}_{i = 1}^{m} $ where $ {X}_{i} \in \mathbb{R}^{p \times q} $. Given a Linear Operator $ \mathcal{A}: \mathbb{R}^{p \times q} \to \mathbb{R}^{k \times l} $ and we have $ {B}_{i} = \mathcal{A} \left( {X}_{i} \right) $.

In order to create the matrix form we should represent the Linear Operator by the matrix $ A \in \mathbb{R}^{\left( k l \right) \times \left( p q \right)} $ which is the matrix which applies the operator on the Column Stacked (See Vectorization Operator) vectors. So we have $ \boldsymbol{x}_{i} = \operatorname{Vec} \left( {X}_{i} \right) $ and $ \boldsymbol{b}_{i} = \operatorname{Vec} \left( {B}_{i} \right) $.

By defining $ X = \left[ \boldsymbol{x}_{1}, \boldsymbol{x}_{2}, \ldots \boldsymbol{x}_{m} \right] $ and $ B = \left[ \boldsymbol{b}_{1}, \boldsymbol{b}_{2}, \ldots \boldsymbol{b}_{m} \right] $ we have $ A X = B $.

With this definition everything will work as expected. Same logic will work for $ F $.

Remark

If you want $ F $ to work along the rows of $ X $ then set $ Z = F {X}^{T} $.

ADMM for the Vector Case

First, pay attention that the ADMM works on vectors (You can work with Matrices but then you need to update the Prox operations accordingly.

So $ \boldsymbol{x} \in \mathbb{R}^{n}, A \in \mathbb{R}^{m \times n}, \boldsymbol{b} \in \mathbb{R}^{m} $ and $ F \in \mathbb{R}^{o \times n} $ where $ o $ is set to match the operation of TV on the image columns and rows which $ \boldsymbol{x} $ is a column stack of.

So now the terms do align propely: $ {\left( {A}^{T} A + p {F}^{T} F \right)}^{-1} \in \mathbb{R}^{n \times n} $ and $ \left( {A}^{T} \boldsymbol{b} + p {F}^{T} \left( \boldsymbol{x} - \boldsymbol{u} \right) \right) \in \mathbb{R}^{n} $. So the calculation of $ \boldsymbol{x}^{k} $ is well defined.

Given that $ \boldsymbol{z} \in \mathbb{R}^{o} $ and $ \boldsymbol{u} \in \mathbb{R}^{o} $ their calculation is also well defined.

I guess something in your code doesn't match. But if the code have the terms $ A \boldsymbol{x} - \boldsymbol{b} $ and $ F \boldsymbol{x} $ well defined then everything else will work.


Conversion between DCT to DFT frequency domains It is well known that the DFT can be interpreted as a filter bank. The same can be done for the DCT, for instance Fig. 3.29, pag. 129 in [1] shows a plot of the frequency response of a DCT filter bank. I am wondering how this plot (and similar plots for other transforms like DST, etc.) was obtained. For the DFT, we have the prototype filter $H_0(z)=\sum_{n=0}^{N-1}z^{-n}$. The frequency response of the $k$-th DFT bin is $H_k(z)=H_0(z\, e^{-j\frac{2\pi}{N}k})$. Using the transformation $z=e^{j\omega}$ we can readily find the frequency response for all bins (see the plot below).

As far as I know there is no z-transform equivalent for the DCT, so it's not clear to me how that plot was obtained. Any ideas?

[1] Malvar, H.S., Signal processing with lapped transforms, 1992, Artech House

Both are linear transforms. Both can be expressed as a matrix. These matrices  are invertible. So you can have the product of one matrix with the inverse of the other. In practice, you could not actually calculate inverse matrices, so you'd have to stick with solvers. However, such products are a little op-art. Coursey GNU octave.