A few questions about Slepian and generalized gaussian windows I'm trying to add documentation for all the window functions in scipy.signal, and I'm stuck on the Slepian (same as DPSS?) and Generalized Gaussian windows, which I'd never heard of before.

There are two variables that are shape parameters of some type, p in the generalized Gaussian, and width in the Slepian.  (sig appears to be sigma, the standard deviation.)

2 questions:

Instead of me reverse-engineering and guessing, can anyone explain what these variables are called and what they do?
Can you explain what these windows are useful for or where they are used?

def general_gaussian(M, p, sig, sym=True):
    """Return a window with a generalized Gaussian shape.

    The Gaussian shape is defined as ``exp(-0.5*(x/sig)**(2*p))``, the
    half-power point is at ``(2*log(2)))**(1/(2*p)) * sig``.

    """
    if M < 1:
        return np.array([])
    if M == 1:
        return np.ones(1, 'd')
    odd = M % 2
    if not sym and not odd:
        M = M + 1
    n = np.arange(0, M) - (M - 1.0) / 2.0
    w = np.exp(-0.5 * (n / sig) ** (2 * p))
    if not sym and not odd:
        w = w[:-1]
    return w

def slepian(M, width, sym=True):
    """Return the M-point slepian window.

    """
    if (M * width > 27.38):
        raise ValueError("Cannot reliably obtain slepian sequences for"
              " M*width > 27.38.")
    if M < 1:
        return np.array([])
    if M == 1:
        return np.ones(1, 'd')
    odd = M % 2
    if not sym and not odd:
        M = M + 1

    twoF = width / 2.0
    alpha = (M - 1) / 2.0
    m = np.arange(0, M) - alpha
    n = m[:, np.newaxis]
    k = m[np.newaxis, :]
    AF = twoF * special.sinc(twoF * (n - k))
    [lam, vec] = linalg.eig(AF)
    ind = np.argmax(abs(lam), axis=-1)
    w = np.abs(vec[:, ind])
    w = w / max(w)

    if not sym and not odd:
        w = w[:-1]
    return w

Possible matches:

nipy's dpss_windows function uses NW, the "standardized half bandwidth corresponding to 2NW = BW*f0 = BW*N/dt but with dt taken as 1"

Matlab's dpss uses time_halfbandwidth Is this the same window?  Is time_halfbandwidth the same thing as width?

This DPSS definition has $ \omega_c$ "the desired main-lobe cut-off frequency in radians per second".

Generalized normal distribution has β (equal to twice p?) which is just called a shape parameter, with normal distribution for β = 1 and Laplace distribution for β = 2.

The Slepian sequences are a family of functions.  Most algorithms calculate 2*NW - 1  sequences at once for a given NW.  N is the number of points in the sequence and W decides half the width of the mainlobe in the frequency domain for the Fourier transform of a given Slepian sequence.  Typically you would use an NW of 3 or 4 for your signal processing.

In scipy, they are asking for the $N$ (m in the Python code) and the $W$ (width in the Python code) parameter separately while in matlab you input the time-bandwidth product $NW$ as a single parameter.  This makes sense because you're typically calculating the Slepian sequences for a fixed window size, $N$.

If you are estimating the power spectrum of a stationary time series the DPSS are the set of windows you should use.

The generalized Gaussian function returns a Gaussian-like function raised to successively higher powers depending on the p parameter.  As p is raised to successively higher powers the generalized Gaussian becomes narrower in the time domain.  The nice property of a Gaussian is that it is it's own Fourier transform and it is the function that achieves the limit concerning the uncertainty principle.  A Gaussian function may be useful if you would like to compute a short-time Fourier transform or spectrogram as an estimator for a time-varying power spectrum of a non-stationary time series.


Meaning of Phase response of a filter? In simple words? What is meant by phase response of a filter?  Often it is associated with restoration of signal shape,so if a filter has linear phase response, its output looks like input otherwise output is distorted if the filter has non linear phase response

Example of linear phase response filter is FIR filter and IIR filter is example of non linear phase response

But what is meant linear or non linear phase response? Please kindly explain in simple words preferably with an example containing figures or graphs for demonstration

Here is the same experiment done in my blog: https://poweidsplearningpath.blogspot.com/2020/04/chapter-51-meaning-of-general-linear.html

First, there is a tiny mistake in the question. Not all FIR filters have the properties of Linear Phase. Only the four types of FIR possess the properties. 1 By contrast, all of the IIR filters are not linear phases.

For me, one better description of general linear phase (GLP) is constant group delay. By definition, group delay is negative derivative of phase 2 (In fact, the detail phase/delay relationship can be derived but I suggest we just accept them.). Then, Derivative (group delay) of a linear (phase) is constant and vise versa.

Then, let's move on to the meaning of group delay.
group delay of a frequency represents the delay unit of the filter to that frequency.
So, the filter may treat different frequencies with different delay unit. For an extreme bad example of an non linear filter, the input signal 'do re mi' may becomes 're mi do' in the output. A GLP filter can guarantee that such wield condition will never happen.

Here I wrote an example. The example comes from the Chapter 5.1.2 in the bible of DSP3 and I just implemented the sample.

First, a IIR Filter with phase response like this is given.

Fig. 1.

Here is the group delay (negative derivative of phase response) and the magnitude response. Please note that I delay the frequency in +-0.2 pi for about 150 units. BTW, the filter is a low pass filter so signal higher than 0.8pi is expected to be filtered.

Fig. 2.

Then, let's input a test signal like 'do re mi'. The signal x[n] are '0.8pi, 0.2pi, 0.4pi' in order. The corresponding frequency response is also provided.

Fig. 3.

And here is the output signal. The signal becomes 'empty, 0.4pi, 0.2pi'. The signal component with 0.8pi is filtered out as expected.

Fig. 4.

To make thing more clear, here I points the number of Fig. 3 and Fig. 4 together. For 0.2pi component, the group delay is about 6.39 units but the group delay of 0.2pi component is about 153 unit. The output signal can confirm the prediction from group delay response. That is why the 0.2pi component becomes that last in the output.

Fig. 5.

In summary,

linear phase equals to constant group delay.
GLP FIR filter can guarantee such scenario will never happen. But IIR can never achieve GLP. (But with the same mag frequency spectrm requirement, IIR usually can achieves the spec with lower delay (but not constant) in comparison with FIR.)

Reference:

FIR filter with linear phase, 4 types
https://en.wikipedia.org/wiki/Group_delay_and_phase_delay
A. Oppenheim and R. Schafer, Discrete Time Signal Processing 3rd. 2009

Matlab code

%% System
% H1[z]
b1 = conv([1 -.98*exp(j*.8*pi)],[1 -.98*exp(-j*.8*pi)]);
a1 = conv([1 -.8*exp(j*.4*pi)],[1 -.8*exp(-j*.4*pi)]);
H1 = tf(b1,a1,-1,'Variable','z^-1');

% H2[z]
H2 = tf(1,1,-1,'Variable','z^-1');
for k = 1:4
    ck = 0.95*exp(j*(0.15*pi+0.02*pi*k));
    ck_conj = conj(ck);
    b_tmp = conv([ck_conj -1],[ck -1]);
    b_tmp = conv(b_tmp,b_tmp);
    a_tmp = conv([1 -1*ck],[1 -1*ck_conj]);
    a_tmp = conv(a_tmp,a_tmp);
    H_tmp = tf(b_tmp,a_tmp,-1,'Variable','z^-1');

    H2 = series(H2,H_tmp);
end

% H[z]
H = series(H1,H2);

% Zero-Pole Plot, Fig. 5.2
[b_h,a_h] = tfdata(H );
b_h = cell2mat(b_h);
a_h = cell2mat(a_h);

figure;
zplane(b_h,a_h);
suptitle('Zero-Pole Plot, Fig 5.2');

% System Response.
L=1000;
dw=2*pi/L;
w = -pi:dw:pi-dw;
HH=freqz(b_h,a_h,w);
mag=abs(HH);
phase=angle(HH);

% Fig. 5.3
figure;
subplot(2,1,1);
plot(w,phase);
xticks([-pi -0.8*pi -0.6*pi -0.4*pi -0.2*pi 0 0.2*pi 0.4*pi 0.6*pi 0.8*pi pi]);
xticklabels({'-\pi','-0.8\pi','-0.6\pi','-0.4\pi','-0.2\pi','0','0.2\pi','0.4\pi','0.6\pi','0.8\pi','\pi'});
xlim([-pi pi]);
yticks([-4 -2 -0 2 4]);
ylabel('ARG[H(e^(^j^w^)]');
xlabel('w');
title('Phase response');

subplot(2,1,2);
plot(w,unwrap(phase));
xticks([-pi -0.8*pi -0.6*pi -0.4*pi -0.2*pi 0 0.2*pi 0.4*pi 0.6*pi 0.8*pi pi]);
xticklabels({'-\pi','-0.8\pi','-0.6\pi','-0.4\pi','-0.2\pi','0','0.2\pi','0.4\pi','0.6\pi','0.8\pi','\pi'});
xlim([-pi pi]);
ylabel('arg[H(e^(^j^w^)]');
xlabel('w');
title('Unwrap Phase response');
suptitle('ARG/arg Plot, Fig 5.3');

% Fig. 5.4
figure;
subplot(2,1,1);
plot(w(1:end-1),-1*diff(unwrap(phase))./diff(w));
xticks([-pi -0.8*pi -0.6*pi -0.4*pi -0.2*pi 0 0.2*pi 0.4*pi 0.6*pi 0.8*pi pi]);
xticklabels({'-\pi','-0.8\pi','-0.6\pi','-0.4\pi','-0.2\pi','0','0.2\pi','0.4\pi','0.6\pi','0.8\pi','\pi'});
xlim([-pi pi]);
ylabel('grd[H(e^(^j^w^)]');
title('Group Delay');

subplot(2,1,2);
plot(w,mag);
xticks([-pi -0.8*pi -0.6*pi -0.4*pi -0.2*pi 0 0.2*pi 0.4*pi 0.6*pi 0.8*pi pi]);
xticklabels({'-\pi','-0.8\pi','-0.6\pi','-0.4\pi','-0.2\pi','0','0.2\pi','0.4\pi','0.6\pi','0.8\pi','\pi'});
xlim([-pi pi]);
ylabel('|H(e^(^j^w^)|');
title('Magnitude response');
suptitle('GD/mag Plot, Fig 5.4');

%% Signal
M = 60;
n = 0:M;
w = 0.54-0.46*cos(2*pi*n/M);

N = 512;
x1 = zeros(1,N);
x2 = zeros(1,N);
x3 = zeros(1,N);
dw = 2*pi/N;
w_freq = -pi:dw:pi-dw;

for i = 0:M

    x1(i+M) = w(i+1)*cos(0.2*pi*i);
    x2(i+2*M-1) = w(i+1)*cos(0.4*pi*i-pi/2);
    x3(i+1) = w(i+1)*cos(0.8*pi*i+pi/5);
end
x = x1+x2+x3;
X = abs(fft(x));
X = fftshift(X);

% Fig. 5.5
figure;
subplot(2,1,1);
plot(x);
title('x[n]');
xlim([0,300]);
subplot(2,1,2);
plot(w_freq,X);
xticks([-pi -0.8*pi -0.6*pi -0.4*pi -0.2*pi 0 0.2*pi 0.4*pi 0.6*pi 0.8*pi pi]);
xticklabels({'-\pi','-0.8\pi','-0.6\pi','-0.4\pi','-0.2\pi','0','0.2\pi','0.4\pi','0.6\pi','0.8\pi','\pi'});
xlim([-pi pi]);
ylabel('|H(e^(^j^w^)|');
title('DTFT of X');
suptitle('Input time/Freq., Fig 5.5');

%% Output
y = filter(b_h,a_h,x);

% Fig. 5.6
figure;
plot(y);
xlim([0,300]);
xlabel('n');
title('output y[n], Fig 5.6');

%= Compre the Delay sample point.
figure;
subplot(2,1,1);
plot(x);
xlim([0,300]);
xlabel('n');
ylabel('x[n]');
title('input');
subplot(2,1,2);
plot(y);
xlim([0,300]);
xlabel('n');
ylabel('y[n]');
title('output');


Bounds of the derivative of a bounded band-limited function Let $f(t)$ be a function with properties:

$$\begin{array}{ll}
t\in\mathbf{R}&t\text{ is in reals}\\
f(t)\in\mathbf{R}\text{ for all } t&f(t)\text{ is in reals}\\
|f(t)|<A\text{ for all }t&\text{absolute value of }f(t)\text{ is bounded above by }A\\
\int_{-\infty}^{\infty} f(t) \ e^{- i \omega t} \ {\rm d}t = 0\text{ for all }|\omega|\ge B&f(t)\text{ is band-limited by frequency B in radians}
\end{array}$$

Given $A$ and $B,$ what is the tight upper bound for $|f'(t)|,$ the absolute value of the derivative of the function?

Nothing else shall be assumed about $f(t)$ than what has been stated above. The bound should accommodate for this uncertainty.

For a sinusoid of amplitude $A$ and frequency $B,$ the maximum absolute value of the derivative is $AB.$ I wonder if this is an upper bound, and in that case also the tight upper bound. Or maybe a non-sinusoidal function has a steeper slope.

You'll be interested in Bernstein's inequality, which I first learned about in Lapidoth, A Foundation in Digital Communication (page 92).

With a well-behaved signal $f(t)$ as you defined it above (in particular, $f(t)$ is integrable and bandlimited to $B\,\text{Hz}$, and $\text{sup}\,|f(t)| = A$), then $$\left|\frac{\text{d}f(t)}{\text{d}t}\right| \leq 2AB\pi.
$$

Note that the original result by Bernstein established a bound of $4AB\pi$; later, that bound was tightened to $2AB\pi$.

I have spent some time reading Zygmund's "Trigonometric Series"; all I'll say is that it is the perfect remedy for those under the impression that they know trigonometry. A full understanding of the proof is beyond my mathematical skill, but I think I can highlight the main points.

First, what Zygmund calls Bernstein's inequality is a more limited result. Given the trigonometric polynomial $$T(x) = \sum_{-\infty}^\infty c_k e^{jkx}$$ (with real $x$), then $$\max_x |T'(x)| \leq n \max_x |T(x)|$$ with strict inequality unless $T$ is a monomial $A \cos(nx+\alpha)$.

To generalize this we need a preliminay result. Consider a function $F$ that is in $\text{E}^\pi$ and in $\text{L}^2$. ($\text{E}^\sigma$ is the class of integral functions of type at most $\sigma$ -- this is one of the places where my math starts to fray at the edges. My understanding is that this is a mathematically rigorous way of stating that $f=\text{IFT}\lbrace F \rbrace$ has bandwidth $\sigma$.)

For any such $F$ we have the interpolation formula $$F(z) = \frac{\sin(\pi z)}{\pi}F_1(z),$$ where $z$ is complex and $$F_1(z) = F'(0) + \frac{F(0)}{\pi} + \sum_{n=-\infty}^\infty {^\prime} (-1)^nF(n) \left( \frac{1}{z-n}+\frac{1}{n} \right).$$ (This is theorem 7.19.)

Now we can state the main theorem. If:

$F$ is in $\text{E}^\sigma$ with $\sigma>0$
$F$ is bounded on the real axis
$M=\sup |F(x)|$ for real $x$

then $$|F'(x)| \leq \sigma M$$ with equality possible iff $F(z) = a e^{j\sigma z} + b e{-j\sigma x}$ for arbitrary $a,b$. We suppose that $\sigma=\pi$ (otherwise we take $F(z\pi/\sigma)$ instead of $F(z)$.)

To prove this, we write the derivative of $F$ using the interpolation formula above: $$F'(x) = F_1(x)\cos(\pi x)+\frac{\sin(\pi x)}{\pi} \sum_{n=-\infty}^\infty \frac{(-1)^nF(n)}{(x-n)^2}.$$ Setting $x=1/2$ we get $$F'(1/2) = \frac{4}{\pi} \sum_{n=-\infty}^\infty \frac{(-1)^nF(n)}{(2n-1)^2}$$ which implies $$|F'(1/2)| \leq \frac{4}{\pi} \sum_{n=-\infty}^\infty \frac{1}{(2n-1)^2} = \frac{4M\pi^2}{4\pi} = M\pi.$$

Now we need a nice little trick: Take an arbitrary $x_0$ and define $G(z) = F(x_0+z-1/2)$. Then, $$|F'(x_0)| = |G'(1/2)| \leq M\pi.$$

(TODO: Show the proof for the case of equality. Define $\sum \prime$.)


Is there a way to obtain the original signal (stationary process) from its combination through filtering (matlab) and crosscorrelation? I have a stationary process $w_1(t)$, white in band $B=[-2, 2] KHz$, and another process: $x(t)=w_1(t)-w_1(t+t_0)$, where $t_0=250\mu s$.

I want to re-obtain $w_1(t)$ by filtering $x(t)$ through $h(t)$: $w_1(t)=x(t)*h(t)$ (where $*$ is the convolution operator).

I know that $r_{w_1x}=r_x*h$, ($r_{w_1x}$ is cross-correlation and $r_x$ is autocorrelation)

from this one I can conclude that $S_{yx}(f)=S_x(f)\cdot H(f)$, so my filter will be $H(f)=\frac{S_{yx}(f)}{S_x(f)}$. ($S$ is the Power Spectral Density)

Is this right? Then from antitransforming I should obtain the time domain representation of my filter.

My problem occurs when I try to calculate it with Matlab:

I have already the autocorrelation of x, and its Power Spectral Density, but I can't figure out how to calculate the cross-correlation of my processes. For example $w_1=randn(8000,1000);$, then filtered, but it's still 1000 representations of a white gaussian noise, and so is $x$.

I thought about doing the mean of all 1000 crosscorrelations and then transform with fft:

Rwx=0;
for i=1:1000
    Rwx=Rwx+xcorr(w1(:,i),x(:,i));
end
Rwx=Rwx/1000;
Swx=fftshift(fft(Rwx));

But this doesn't seem to work.

Do you have an idea on how to solve this?

I think that (at least in general) (and if I understand your question correctly) this is not possible. As an example, assume that $w_1(t)$ contains a harmonic component of period $t_0$. The signal $x(t)$ will not contain any trace of this component. As a result no amount of filtering $x(t)$ will be able to recover the original signal.


Signal matching: How to quantify the difference between two signals of varying frequency I have two waves, and I would like to quantify the difference between them.

They are exponential sine waves where the frequency decreases with time.

I'm interested in the difference in decay constant between the signals. I thought maybe I could use some kind of FFT method to get the envelope, but I don't know what to do with the varying frequency.

For a relatively easy answer you can use a short-time Fourier transform and use the largest value as the estimate of your sinusoid frequency.  The basic idea is that you do short Fourier transforms at different points in time to see how it changes.  The problem with that is that you have a tradeoff between time resolution and frequency resolution.  Shortening the Fourier transform improves the time resolution at the expense of the frequency resolution.  Conversely, making the Fourier transform longer improves the frequency resolution at the expense of the time resolution.

The frequency resolution can be improved by using quadratic interpolation to get a finer estimate of the true sinusoid peak.

Another way to go that can get you very good frequency estimation at every point in time is with a phase-locked loop.  Phase-locked loops are a little trickier to set up though, because you have to get the feedback gains right.  There is also a tradeoff in regards to responsiveness (akin to time resolution) and resistance to noise.


which features fundamental frequency is correlated? I'm unable to compute $f_0$ (fundamental frequency) in librosa feature extraction. By ready much in github issues check the second comment. I see that $f_0$ is correlated with some other features (like zero_crosssing) easy to extract from librosa.

What are these features except for that one I give? Is this correlation super high so that they can replace $f_0$?

data sample

site for data

Here is the answer as responded here in github.


Comparing with FIR, why does IIR always gives a bias error, when calculating signal DC? I am trying the DC blocking function in matlab. The algorithm of DC blocking is to find DC of the signal first and then remove the DC value from the signals. Here are the codes I tried, which are sample codes from Matlab:

  t         = (0:0.001:10000)';
  x         = sin(30*pi*t) + 1;

  hDCBlock1 = dsp.DCBlocker('Order', 4);  % IIR is default
  hDCBlock2 = dsp.DCBlocker('Algorithm', 'FIR', 'Length', 100);
  hDCBlock3 = dsp.DCBlocker('Algorithm', 'Subtract mean');
  for idx = 1 : 10
    range = (1:1000) + 1000*(idx-1);
    y1 = step(hDCBlock1, x(range));
    y2 = step(hDCBlock2, x(range));
    y3 = step(hDCBlock3, x(range));
  end
  plot(t(1:1000),x(1:1000),...
    t(1:1000),y1, ...
    t(1:1000),y2, ...
    t(1:1000),y3);
  legend(sprintf('Input DC:%.3f',    mean(x)), ...
    sprintf('IIR DC:%.3f',           mean(y1)), ...
    sprintf('FIR DC:%.3f',           mean(y2)), ...
    sprintf('Subtract mean DC:%.3f', mean(y3)));

I tried to adjust the length of data and order of the IIR, the DC given by IIR always returns a bias error. However the FIR and mean does not, or gives very little errors.

Could anyone explain on this?

i think this is about a well-known limit cycle problem with IIR filters, particularly implemented using fixed-point arithmetic.  this is something i wrote a long time ago. from which approximately half of this article is derived from.

maybe later i will copy my thing over here and dress up the equations nicely with $\LaTeX$.


Equal probability quantizer in Python I am trying to apply Haralick textures to a SAR image (float32). As far as I know, the image first needs to be quantized to a reasonable bit depth prior to calculate the co-ocurrence matrix. In the original paper, Haralick proposes using a "equal probability quantizer". As far as I know is the algorithm employed by Mathlab "imquantize". I have been looking for a pre-existing python/numpy implementation (no success).
I also have reading pre-existing java code from ESA SNAP. As far I know histogram calculations would be straightforward with numpy, but if I understand correctly, every pixel is assigned  a "level/histogram bin" using bisection search, so I do not see any straightforward way to properly vectorize the code.
EDITED after this point ===================
I have made simple implementation in Python + Numpy based on the pre-existing java code from ESA SNAP, and the kind advice of Marcus Müller.
import numpy as np

def eqProbQuant(image, levels=32):
    # Sort the pixels by value
    sorted_image = np.sort(image)
    # Get the pixel count
    pixel_count = sorted_image.shape[0]
    # Get the number of pixels per bin
    samples_per_bin = int(pixel_count/levels)
    # Get locations where the bin would change
    edge_samples = np.arange(levels+1) * samples_per_bin
    # Get the values at those locations (bin edges)
    bin_edges = sorted_image[edge_samples]
    # Use the values to apply quantization
    quantized = np.digitize(image, bin_edges)

    return quantized

If possible, I would like to know if this code would apply the algorithm correctly (I have very little experience with signal processing)

There is any way to assign a a "level/histogram bin" to each pixel without iterating through the entire image?

(emphasis mine)
No, because what you describe is an operation that needs to be applied to each pixel.
Also, a histogram needs to be calculated taking all pixels into account (that's the definition of a histogram!).
Anyway, that's not a problem at all – in SAR processing, this operation is definitely of "a lesser concern" when it comes to computational intensity.

Calculate a list of pixel values sorted by intensity
assign value ranges based on that

Sort the pixels
That's inherently an operation that can't be fully parallelized – in the end, there is one list for the whole image. However, you can parallelize it for the most part: Let your multiple processors (GPU threads, …) sort disjunct parts of the image, then merge the sorted lists in order (and that can be done in cascaded ways, too). Pseudocode:
N := numbers of pixels
I := image considered as linear array of pixels
O := output list of sorted pixel values, same length as I
T := temporary list of list of pixels
P := numbers of processing cores

FOR p IN P PARALLEL:
   T[p][:] := sort(I[p*(N/P) … (p+1)*(N/P)])

merge(T[0…(P-1)])

Quantize
Then, just take the sorted list of values, divide into the desired number of bins of equal size, note down the "starting" value of each bin (i.e. the values at O[0], O[N/n_bins], O[2*N/n_bins]...).
Then, (again, in parallel if desired – there's no data interdependence between pixels) move through your image and replace the pixel values by the bin number.


Getting HRTF from 3-D scan + acoustic simulation This question follows from Modelling propagation of sound wave by particle simulation

This last fortnight I have started experimenting with Binaural audio; it really blows everything else out of the water.

The main reason it hasn't taken off, I think, is that every individual needs to get their own HRTF calculated, which consists of sitting in a soundproof room with microphones in your ear canals, while something makes popping noises at thousands of locations around you.

I'm only aware of one service provider that may calculate your HRTF: http://www.physiol.usyd.edu.au/~simonc/hrtf_rec.htm  and they are in Australia!

I'm wondering whether a more practical method may emerge, which would involve taking a 3-D scan of one's head (maybe using something like http://www.david-3d.com/) and shipping it off for a heavy dose of distributed computing which would return HRTFs for that individual.

I wonder if we may see a day where people have their own HRTF data stored in the cloud, and can enjoy binaural sound on their iPods.

How far away is such a technology?  And other any other contenders for measuring HRTF?

π

PS could we possibly have some more tags like "binaural", "hrtf/hrir"

There are various school of thoughts about this. A German guy named Klaus Genuit came up with a parametric model of HRTFs that's based on the actual physcial geometry of head, shoulder, pinna, etc. He actually now runs a company that craeates binaural products and measurement gear, see for example
http://www.head-acoustics.de/downloads/publications/binaural_technology/DAGA_99_Physical_Modeling_of_individual_HRTFs.pdf

Based on this work you could potentially create individual HRTFs using a some 3D scanning scheme. I do vaguely remember a company coming doing exactly that for a bunch of people at my employer at the time. I didn't think the HRTFs that came out that way were particularly good but I didn't do a lot of work with them either.

Individualized HRTFs only help to a certain extent. Let's say you want to record a concert. You can do this binaurally using a dummy head. However, there is currently no way to convert this recording to your individual HRTFs after the fact, even if you have them.


Simplification of a convolution-like expression I am trying to further reduce the following expression (if there's a nicer way to write it):

$\sum_{k=-\infty}^{+\infty}\sum_{p=-\infty}^{+\infty} \sum_{l=-\infty}^{+\infty} w(k) \varphi_\mathrm{y}(l-p)\varphi_\mathrm{x}(k-n-l+p)$

where $\varphi_\mathrm{y}$ and $\varphi_\mathrm{x}$ are the autocorrelations of the signals $y(i)$ and $x(i)$, respectively.

I know I can express $\sum_{k=-\infty}^{+\infty} w(k) \varphi_\mathrm{x}(k-n-l+p)$ as a convolution (or a correlation, using the evenness of $\varphi_\mathrm{x}$) but I left the expression expanded in case it can help.

If necessary, I could take the simplified case of finite (M-samples-)length signals.

Regards

As you say, since $\varphi_x$ is an autocorrelation we have that
$$\sum_{k=-\infty}^{\infty} w(k) \varphi_\mathrm{x}(k-n-l+p) =
(w \star \varphi)\bigr |_{n+l-p} = g(n+l-p)$$
where $g = w \star \varphi$.
So, since $\varphi_y$ is also an autocorrelation, we have
$$\begin{align}
\sum_{l=-\infty}^{\infty} g(n+l-p)\varphi_y(l-p)
&= \sum_{l=-\infty}^{\infty} g(l+n-p)\varphi_y(p-l)\\
&=\sum_{m=-\infty}^{\infty} g(m)\varphi_y(n-m) &
\text{upon setting}~l+n-p=m,\\
&= g \star \varphi_y\bigr|_n = h(n)
\end{align}$$
where $h = g\star\varphi_y = w\star\varphi_x \star\varphi_y$.
So now you are left with $\displaystyle \sum_{p=-\infty}^\infty h(n)$ which
is unbounded. Are you sure that this last sum is with respect to $p$
and not with respect to $n$?


Dominating frequency in a "clear" 50 Hz signal I am trying to find out the dominating frequency of a signal with a frequency of 50 Hz (sampled at 200 Hz - every 5 milliseconds).

The python code I am using to do this is the following (based on this):

import numpy
from numpy import sin
from math import pi

t = numpy.linspace(0, 1, 201)  # 200 Hz sampling rate
y = sin(2*pi*t*50)

fourier = numpy.fft.fft(y)
frequencies = numpy.fft.fftfreq(len(t), 0.005)  # where 0.005 is the inter-sample time difference
positive_frequencies = frequencies[numpy.where(frequencies > 0)]
magnitudes = abs(fourier[numpy.where(frequencies > 0)])  # magnitude spectrum

peak_frequency = numpy.argmax(magnitudes)

The peak frequency I get is 49 Hz, not 50 Hz (what I think I should obtain).

Is this the correct way to calculate the peak frequency? If yes, then can you please explain why this happens? Where is my mistake?

EDIT
Code modified according to SergV's answer:

import numpy
from numpy import sin
from math import pi

Fs=200.
F=50.
t = [i*1./Fs for i in range(200)]
y = sin(2*pi*numpy.array(t)*F)

fourier = numpy.fft.fft(y)
frequencies = numpy.fft.fftfreq(len(t), 0.005)  # where 0.005 is the inter-sample time difference
positive_frequencies = frequencies[numpy.where(frequencies > 0)]
magnitudes = abs(fourier[numpy.where(frequencies > 0)])  # magnitude spectrum

peak_frequency = numpy.argmax(magnitudes)

The calculated peak frequency is still 49 Hz...

EDIT 2:
It seems that the problem was that the positive frequencies were not including 0. After changing the code to

positive_frequencies = frequencies[numpy.where(frequencies >= 0)]
magnitudes = abs(fourier[numpy.where(frequencies >= 0)])

the calculated dominating frequency was 50 Hz (as expected).

This is the magnitude spectrum I get: http://imgur.com/JxhnNc5 It shows that the dominating frequency is 50 Hz. Interestingly, even when, because of the bug, the calculated peak frequency was 49 Hz, the magnitude spectrum was showing the correct dominating frequency (50 Hz)

The size of $y$ is 201. So frequency bin of fft is not 1 Hz. Change code for calculating $t$:

Fs=200.
F=50.
t = [i*1./Fs for i in range(200)]

Now you have only 1 peak on 50 Hz.

P.S.

Plot magnitude spectrum
Try to modify F (49, 49.5, 49.67). Enjoy!

Welcome to DSP and new questions!


How can one put information within a signal waveform and make it to traverse through different kind of medium? This is a soft question. There are different types of signals
like analog ,discrete or digital one.

But my question is how can one put information within a signal waveform and make it to traverse through different kind of physical medium, electronic circuits,devices etc?

Moreover ,in physics,we have Newtons laws for conservation of mass and energy . Now,is there any similar law in case of signals?

Maybe what you want is information theory. Information can be quantified by Shannon entropy. An analog signal channel always has some level of noise which limits the rate of information that can be sent through it. That limit is called the Shannon limit, or more specifically for a given channel bandwidth the Shannon–Hartley theorem.


Why is the frequency measured in radians in filter design? This is my first time posting a question, so sorry if my articulation is not proper.

I'v just started out with FIR filter design and I see the domain being -$\pi$ to $\pi$. Now I understand that discrete frequencies repeat after an interval of 2$\pi$ but consider this. I need to design a low pass FIR filter with a cutoff frequency of 100$Hz$. How can I specify this in terms of radians? Is it using $\omega$=$2$$\pi$$f$? But that would be outside the domain. As you can see, I am very confused about this concept so any help would be appreciated.

The value of $2\pi$ radians corresponds to the sampling frequency. So the normalized frequency in radians is given by

$$\omega=\frac{2\pi f}{f_s}\tag{1}$$

where $f$ is the actual frequency in Hertz, and $f_s$ is the sampling frequency in Hertz.


What base image size to choose to assist easier downsizing? Android applications need image at different resolutions to support different screen densities and pixel size.

I am planning to provide only one image with high pixel size and let the underlying platform do the resizing. When we shoot an image with a professional camera, let's say I get image size of $5000\times 4000$ pixels. Which is obviously too high for mobile devices. Now to reduce the size of the file (to be bundelled with namtive mobile application) I have two option:

Keep the original image pixel size and perform high compresson at jpeg level
Downsize the image to more reasonable mobile screen upper limits say $2000\times 1600$ and keep jpeg compression low?

Which option is better quality wise considering target file size to be same?

There is no objectively correct answer, there are too many variables to consider such as the nature of the image you are compressing (how much entropy/how 'compressible'?), how much you would need to reduce the jpeg quality to achieve the desired filesize etc... If it is possible to use a lossless compression format then I would consider that as a solution.

If you want to take an empirical approach you can take some representative images, go through the process of converting them using your two methods, then compare the results. You can use a tool like compare for a qualitative metric or visual comparison.


Optimal segment length for coherence estimation I have a question about the computation of coherence between two signals with Welch's method.
The two signals are relatively short (e.g. 256 samples) and I would like as correct an estimate of shared frequencies between the signals as possible.
I am aware that I need to average over more than one segment of data for meaningful results and that overlapping segments will improve accuracy. Computational effort in the computation is not an issue, so I am using $noverlap = nseg-1$. Acquiring more samples is not possible.

The coherence is computed as $C_{xy} = \frac{|G_{xy}|^2}{G_{xx} G_{yy}}$ where $G_{xy}$ is the Cross-Spectral Density and $G_{xx}$ and $G_{yy}$ are the Power Spectral Densities of both signals, respectively.

My question is:
What is the optimal length of segments?
nseg=ndata (no averaging) leads to $C_{xy} = 1$ whereas a too small value of nseg leads to insufficient frequency resolution.
Is there any value for nseg for which an optimal compromise between sufficient averaging and frequency resolution can be found?

Both theoretical and best-practice suggestions would be appreciated!

There is good news and bad news. You are actually showing magnitude squared coherence(MSC). There is a very complete treatment of MSC, including closed form equations for the sample MSC as a function of independent samples.  There is also a discussion on overlap in:

    G. C. Carter, “Coherence and time delay estimation,”
 Proceedings of the IEEE, vol. 75, no. 2, pp. 236–255, 1987

The bad news is in the plot, which is the sample pdf for MSC where the true MSC is .3, in a single bin where the legend denotes the number of independent terms. If you want 256 frequency bins,  you really need around $256^2$ data points.

You might want to look up the paper and note any recent citations where there is some scheme to reduce the amount of data you need,  but 256 points is not anywhere near enough.


Where did we get the DC term of the Accumulator from DTFT? Define $y[n]:=\displaystyle\sum_{m=-\infty}^{n}x[m]$. The DTFT is found as follows:
\begin{align*}
y[n]&=\sum_{m=-\infty}^{n}x[m]\\ \\
&=\sum_{m=-\infty}^{n-1}x[m]+x[n]\\ \\
&=y[n-1]+x[n]\\ \\
\implies x[n]&=y[n]-y[n-1]
\end{align*}
Therefore, applying DTFT :
\begin{align*}
\mathcal{F}\{x\}[n]=X(e^{j\omega})&=\mathcal{F}\{y[n]-y[n-1]\}\\ \\
&=(1-e^{-j\omega})Y(e^{j\omega})
\end{align*}
Hence :
$$
Y(e^{j\omega})=\frac{1}{1-e^{-j\omega}}X(e^{j\omega})
$$
However, the correct answer is :
$$
\frac{1}{1-e^{-j\omega}}X(e^{j\omega})+\color{red}{\pi X(e^{j0})\sum_{k=-\infty}^{\infty}\delta(\omega-2k\pi)}
$$
My question is where did the DC term come from and why didn't it work in my derivation?

Note that the equation
$$x[n]=y[n]-y[n-1]\tag{1}$$
doesn't uniquely determine the sequence $y[n]$. If some $y[n]$ satisfies $(1)$, so does $y[n]+c$ with some constant $c$. Consequently, Eq. $(1)$ determines $y[n]$ only up to a constant, which corresponds to a DC term in the frequency domain.
Hence, the equation
$$Y(e^{j\omega})=\frac{X(e^{j\omega})}{1-e^{-j\omega}}\tag{2}$$
is only valid for $\omega\neq 0$, or, more generally, for $\omega\neq 2\pi k$, $k\in\mathbb{Z}$.
The most straightforward way to derive the DTFT of
$$y[n]=\sum_{m=-\infty}^nx[m]\tag{3}$$
is to realize that $y[n]$ is the convolution of $x[n]$ with the unit step sequence $u[n]$, and, consequently, the DTFT of $y[n]$ is given by
$$Y(e^{j\omega})=X(e^{j\omega})U(e^{j\omega})\tag{4}$$
where $U(e^{j\omega})$ is the DTFT of the unit step $u[n]$:
$$U(e^{j\omega})=\pi\delta(\omega)+\frac{1}{1-e^{-j\omega}}\tag{5}$$
Several ways to derive the result $(5)$ are discussed in the answers to this and this question.
Note that in $(5)$ I use $\delta(\omega)$ instead of $\sum_k\delta(\omega-2\pi k)$ to avoid cluttered notation. The fact that the DTFT is always $2\pi$-periodic is understood.


How to discover an image clone I have an image depicting a results of gel electrophoresis (it does not really matter - it just means: no bodies, faces, predictible objects) and I have a strong suspicion that a part of the image background was copied (probably Photoshop Clone Stamp) to cover an unwanted peak. The background is really just noise, but it should be almost perfect copy. I am not entirely sure how big is the copied area so just straightforward cross-correlation with an image part won't work.

How do I prove that? I prefer Matlab for such tasks.

Maybe this is an appropriate paper for you: Exposing Digital Forgeries in Scientific Images, by Hany Farid (freely available, e.g. at http://www.ists.dartmouth.edu/library/199.pdf).


Output of MATLAB's two-dimensional DWT function $\tt dwt2$ This is a conceptual doubt about the two-dimensional DWT. I am trying to understand the decomposition step of a two-dimensional DWT.

In the MATLAB explanation, the two-dimensional DWT generates 4 new images (cA, cH, cV, and cD) that represent the approximation, the horizontal details, the vertical details and the diagonal details.

I am also reading the explanation about this two-dimensional DWT. But the author shows the results as $\text{LL, LH, HL and HH}$.

I want to know if it corresponds to the same thing as cA, cH, cV, and cD. Is the following comparison correct?
\begin{align}
&\texttt{cA} =\text{LL} =\text{approximation}\\
&\texttt{cH} =\text{LH} =\text{horizontal}\\
&\texttt{cV} =\text{HL} =\text{vertical}\\
&\texttt{cD} =\text{HH} =\text{diagonal}
\end{align}

I'm also a little confused about the $\text{LH}$ and $\text{HL}$ results. In some posts, I have read that the $\text{LH}$ represents the horizontal details and in other, I have read that it represents the vertical details.

For example, as answered here , the user explains that:

  In the wavelet terminology, LL is the approxmation image, LH is the vertical details, HL is the horizontal details, and HH is the diagonal details.

But in this article the author says:

  The LL band corresponds roughly to a down-sampled (by a factor of two) version of the original image. The LH band tends to preserve localized horizontal features, while the HL band tends to preserve localized vertical features in the original image. Finally, the HH band tends to isolate localized high-frequency point features in the image.

Can someone clarify this for me?

The confusion is natural as the notations and practice are misleading, and differ across the different works on 2D wavelets.

First, the order of operations ("row then column" or "column then row") is important. In most cases (to the best of my knowledge), LH means low-pass on rows, then high-pass on columns (so, giving horizontal edges). But in some contexts, people use filtering tools on "columnized signals", so the first filter works on columns, then the image matrix is transposed, the second filter pass is applied, and the result is transposed back.
Second, from an math operator point of view, LH(I) is more often understood as applying the high-pass first, then the low-pass, enhancing the confusion.
Third, one may think about the active derivative: a vertical derivation yields the horizontal features.

In the following image, the rightmost letter seems to denote the first filtering operation:

while the opposite convention is used in second image:

This might stems from a cultural and technical background: does one read images left-right then top-bottom, or  top-bottom then left-right?

The best is to test your implementation, look at the resulting images, and use a consistent notation.


Separate static echos from acoustic signal (ICA?) I am using acoustic waves to measure the thickness of material. My problem is that the acoustical lens produces static echos in the signal. I try to separate both because there is a short interference between both.

The frequency of the material echo and the static echo has nearly the same frequency. So spectral subtraction is not a valid method to separate the signals.

I will try the following:
The first measurement is without any object. So I am able to measure the static echo of the acoustic system. The second measurement I collect the signal with the object to measure.

My idea:
So I got two signals from different signal sources. Maybe I can separate both with Independent Component Analysis (ICA).

Do you think that this solution will work? Or is there a better method to separate static and echo signals (with noise) out the measured data.

Thanks for your help.

This probably should be a comment, but I didn't earn that option yet. You can try measuring instrument impulse response. If I understand correctly, your system is essentially one-dimensional. Then if you use thin plate or plate of known thickness instead of thick object/material, your output signal will be an approximation of impulse response. Using it as deconvolution kernel will allow you to measure response specific to material, or object.


Image distortion and alignment on 3D mesh / pointcloud Currently I have a fairly accurate (50 µm) and large dataset from a 3D scanner and a set of high resolution pictures. Sadly the 3D scanner (a Nikon/Metris arm-based CMM that was hanging around the lab) doesn't record colour data, but this would be exactly what I need to extract the actual features I'm interested in. I started out by correcting the lens distortion for the camera, and also accounted for the non-linear distortion. What prevents me from doing this manually in software like Meshlab or Geomagic is that I have an extremely large set of objects to scan, which makes it highly unrealistic to do it by hand.

I took the liberty to add a series of reference points which are visible to both the 3D scanner and conventional camera. I am able to detect these features on both the point cloud and images automatically. Sadly I must admit I underestimated a big part of this task: I now wish to automatically align the image versus the point cloud based on these reference points and project the image onto the point cloud or mesh. The reference points from the point cloud are in a XYZ format, while those from the image are 2D coordinates in the plane. They are already matched, so that presumably - according to a fellow Ph.D. student - difficult part is already done. And while I can think of a few ways to align, distort, and project the images this would most likely be a large investment of time from my part. I was wondering if any pre-existing libraries for this exist or at least if there are papers which present a method to do so. Based on the fact that software like Geomagic is able to do it I'd like to avoid reinventing the wheel. However I must humbly admit I can't seem to find them as I am probably using the wrong terminology.

If you're only left with matching 3D points to their corresponding 2D parts, this is a fairly task and called PnP pose estimation. OpenCV has a lot of options to solve this problem. Here is one tutorial, which you could maybe benefit.


How to prove that this impulse response is stable? This is the impulse response:

Can anyone give a detailed method to prove that it is stable or not ?

If you check graph of h(t) which is underdamped oscillations presrnt  in 1st and 4th quadrant only. And integration of this h(t) from limits (-infinity to + infinity )  Is finite value. therefore,  system is stable


Gaussian window and standard deviation I am trying to figure out how to use the matlab gausswin function which constructs a Gaussian window for $N$ samples with a given standard deviation $\sigma$. The function is defined by

w = gausswin(N, alpha)

where alpha is $1/\sigma$. I would like my Gaussian window to contain $\pm 3 \sigma$ in the Gaussian. So for $N$ samples it seems I should do: $6 \sigma = N$ or $\sigma = N/6$. However, here are the results for $N = 100$ and $\sigma = 100/6 \approx 16.666$:

Clearly this Gaussian does not extend to $\pm 3 \sigma$. Can someone tell me what I'm doing wrong?

My guess it is has something to do with some length scale which isn't being properly specified. It seems to me that the alpha input to the gausswin function is some kind of normalized (unitless) standard deviation. Perhaps I'm not accounting for that.

The manual entry does not say that $\alpha = 1/\sigma$. It says $\alpha$ is proportional to the inverse of $\sigma$.

If you look at the example on that page you'll see

stdev = (N-1)/(2*alpha);

which, for your values says

stdev = (100-1)/(2*16.6667) = 825

you should be scaling the relationship by $N-1$.


Cross-correlation of two signals using matlab I am trying to measure the similarity between two signals and I am using cross-correlation to achieve this.

Following the matlab example, I seem to have accomplished that, however, I do not understand the plot that I am getting. In the matlab example the explanation that is given is: The first subplot indicates that the signal and template 1 are less correlated while the high peak in the second subplot indicates that signal is present in the second template.

My question is how does the peak indicate that the signals have more in common? Shouldn't the two signals be more correlated when the distance between them is lower?

Moreover, when I tried it using my two signals I received the following image:  which looks more than strange compared to the image from the example, as I have multiple peaks and a relatively long straight line.

EDIT

I added two example signals that yield a relatively same result. As suggested from the comments, I tried without upsampling the signal and making them equal in length, however, the results weren't appealing as well. Furthermore, if the peak is, for example, at the 10s, does this indicate that the signals are simillar at the tenth second? What about if the peaks are in the negative time zone (t<0s)?

xcorr function (of Matlab) computes the cross-correlation between two sequences $x[n]$ and $y[n]$ of length $M$ each. If $x$ and $y$ are random processes (jointly WSS) then xcorr returns an estimate $$\hat{r}_{xy}[m] = \frac{1}{M} \sum_{n=0}^{M-1} x[n]y[n+m]^*$$ computed over the two random vectors $x$ and $y$, of the theoretical cross-correlation between the two random processes $X$ and $Y$:
$$r_{xy} [m] = E\{ X[n]Y[n+m]^* \} $$ (Note that the practical estimate computation (the sum's range and the scale) may be modified to make it an unbiased estimator and to deal with out of bound entries...)

When the two sequences belong to deterministic signals, then it's called as the deterministic cross-correlation, essentially the same thing computed from a practical point of view, but its interpretation is made carefully depending on the case.

In signal processing theory and signal detection theory, the operations of correlation, convolution and matched filtering merge at some point. Without further details, one can observe that given two sequences $x$ and $y$, the the computation mechanics and hence the results of the three operations, will be the same (with slight modifications in the orientation of the arguments)

Matched filtering, essentially convolves a test signal with a target signal, with the test signal reversed in time. It can be shown that there'll be a unique peak at the matched filter output, if the target signal involves a copy of the test signal inside. From a theoretical point of view, intuition tells that, if the target and test signals are uncorrelated (and orthogonal, for zero mean processes, as a result) their cross-correlations should produce low values in line with inner products between orthogonal vectors. But when the target signal involves a copy of the test signal at some point, then there will be high correlation between the target and test signals close to those shifts, at the output, that will create the peak.

There are subtle details, however, such as if the test signal is contained more than once, if several test signals overlapp, and if the test signal is not contained purely but through some (inevitable) transformations, which make the comparison more difficult or less reliable unless those transformations are known and rolled back properly.

In most practical scenarios, the target signal will consist of noise  or irrelevant signals plus the test signal. And for those lags which compare the test signal with the noise portions, the output will be arbitrary (in the ideal case of when the noise belongs to a zero mean white random process, the output should tend to zero, as the summation would sequentially add and subtract test signal samples which would produce apprx a zero, when the test signal belongs to a zero mean (ergodic) process). And when the comparison is between target signal and test signal for a lag that corresponds to the place of the test signal in the target signal, the computation should produce a large value, as the sum is essentially computing the energy of the test signal; i.i, sum of $x[n]^2$ values.

Note that when the SNR is low; that the noise portions have large amplitudes compared to the test signal embedded in the target signal, then the above paragraph becomes invalidated and the expected peak is burried within random peaks produced as the result of the irrelevant computations between the test signals and noise contained in the target signal.


Specific Frequency Resolution I have an audio signal

SampleRate Fs: 44100 Hz

TotalSamples: 94144 samples

Duration t: 2.1348 s

The frequency resolution is given by  Fs/N where FS is the input signal's sampling rate and N is the number of FFT points used.

If I want to have 10Hz frequency resolution (bins of 10Hz), I should use 4410 FFT-points (44100/4410 = 10) but fft(signal,N) function of matlab specifies that N (4410 in my case) should greater than the signal length (94144 in my case). How should I proceed?

There are some important clarifications out of this question that are good to point out:

The frequency resolution of a block of data is 1/T where T is the length of the data in time (in seconds). Since the sampling rate and block length are related as

$T = N/f_s$

Where $f_s$ is the sampling rate, it follows that the frequency resolution for a block of N samples will be $1/T = f_s/N$.

Zero padding does not change the frequency resolution, it only interpolates more samples of the Discrete Time Fourier Transform. See these posts:

What happens when N increases in N-point DFT

What proportion of a padded FFT should be actual values

When you do fft(signal,N) in Matlab, N must be larger or equal to the signal length and when larger it will simply append zeros, so identical to zero padding.

So to achieve a 10 Hz resolution bandwidth, you either need to reduce the number of samples, as in fft(sig(1:m)), or reduce the sampling rate, but regardless of sample rate, the length of the block in your fft must be 1/10 = 0.1 seconds long.

Note that windowing if performed will expand the frequency resolution; the finest resolution is achieved with no windowing and is $f_s/N$

Refer to this paper by fred harris which details the resolution bandwidth for various windows: fred harris on the use of windowing


invfreqz() and frequency domain filter specification I am trying to go from a time-domain description of a filter, via frequency domain, generating a filter based on that frequency domain description, then seeing how far from the original I end up.

L = 10-1;
h = zeros(L+1,1);
h(2) = 1;
L2 = 1024-1;
H = fft(h, L2+1);
H1 = H(1:(L2/2+2));
W = 2*pi*(0:(L2+1)/2)/(L2+1);
n = 6;
m = 2;
[b,a] = invfreqz(H1,W,n,m);

ir_in = impz(h);
ir_out = impz(b,a);
[H2,W2] = freqz(b, a);

For a simple input like this, I get a warning about singular or badly scaled matrix, and the solution (as viewed in the filter time domain) is quite different from the original:

b =
    0.0000    1.0000   -1.4219    0.2578   -0.0000    0.0000    0.0000

a =
    1.0000   -1.4219    0.2578

The frequency magnitude and phase response of the new filter looks great, but the numerical properties not so much.

I realize that setting m = 2 is kind of silly in this example, but having an automated and generic process would have been nice.

The warning is due to modeling a unit delay which has wide (constant) frequency response thus will have very large number differences (ill conditioned matrices) in the translation between domains and looking for an IIR solution to the FIR problem (note that the algorithm converges to a trivial but valid solution where the denominator coefficients match the numerator-- if the OP set $m$ to be larger, the algorithm will converge to populate those coefficients as well in the numerator). The recommended approach is to start with FIR only, and only introduce IIR if needed. Regardless the result does match the expected very well in both time and frequency as plotted below after running the OP's code:

Frequency Response

Impulse Response

The above is the IIR equivalent to the OP's FIR target given $m=2$ in the call to invfreqz, which would then converge to an IIR solution. To force an FIR solution, set $m=0$ which results in the following for this case:

b = 0 1 0 0 0 0 0
a = 1

(Note the zeros for b were actually +/-3e-17 or lower)

This is a useful approach for fitting to a precise frequency response curve. For general lowpass, bandpass and multiband filter design, I prefer to use the least squares firls() function with simple targets of 1 and 0 for the passbands and stopbands. However if a more elaborate response is desired (in a curve fitting sense), this is the approach I would take:

Start with a very large $n$ and $m=0$ and run invfreqz() to evaluate the resulting impulse response (b coefficients). Assuming the impulse response converges sufficiently to what can be assumed 0 (finite) at the start and end of the filter (evaluate in dB scale as small values can have a significant impact), reduce $n$ to this length and rerun invfrewqz. Evaluate the frequency response using freqz and increase $n$ if needed to meet target (or reduce $n$ until desired maximum error is reached). The result of this is the best FIR only solution that can be achieved. Given the coefficients of the FIR filter is the impulse response, any filter with an impulse response that converges to zero can be realized as an FIR only filter, but implementing as an IIR can reduce the total number of coefficients needed (with all the cost of IIR implementation such as possible instability, round-off error accumulation, dead-beat responses, etc).

If the number of elements in the FIR is unacceptably long such that an IIR solution is desired, increase $m$ and repeat which will result in a smaller $n$ and an IIR solution. (Make $m$ excessively long to start and evaluate b and a on a log scale to determine the number of significant coefficients). If using an IIR solution, it is recommended to then factor the resulting solution into 2nd order sections to minimize round-off errors that can accumulate significantly in higher order IIR filters. (Resulting in 2nd order biquad filters).


Recommendation for good ressources on digital processing and D-QPSK I've been a C++ developer for 10 years. Quite high level. And I am completely new to Digital Signal Processing.

I am trying to implement a Software Defined Radio. Specifically I am trying to receive FM broadcasts, and more importantly DAB+.

I'm on my way to buy a SDR dongle: Basically a box which links an antenna to a tuner and a ADC. The data is then transferred to the computer through a USB cable. Then it can be delivered to my program thank to a library provided by the SDR dongle vendor.

From there my program will need to demodulate the signal. Well I know from Wikipedia that FM encodes the signal using Frequency modulation. And DAB+ is actually using D-QPSK modulation.

I really don't know how can my program process those modulation in order to get the underlying data.

Do you know any good resources over the internet which could help me?

Any good books for newbies?

If those books contains example code, then that'd be a bonus point. But if the book you want to advise are quite "theoretical", then it's quite fine to list them: the goal of all these is for me to understand the whole machinery of Digital Signal Processing, and I do not want to sacrifice theory for practice.

Wow, I'm honored by Matt L. doing what I'm often doing: Referring people to GNU Radio.

The project actually has a list of recommended literature, but I don't know how well that'd fit you. It's probably still worth looking into.

Then, regarding QPSK: Well, it's one of the basic constellations, and you'd probably be best off reading a textbook intro to digital communications – I like Proakis' Digital Communication, personally, and the third edition (which sadly reduced the formula count) is available cheaply used from libraries.

Regarding DAB+: It uses a differential QPSK, but not directly on the air, but to drive the subcarriers of an OFDM signal, if I'm not mistaken; so you should understand both, as well as (at least superficially) the channel coding used. This again calls for some textbook knowledge!

There's already a Free and Open Source DAB+ transmitter implementation – opendigitalradio is the home of that.

There used to be a GNU Radio DAB (+) transceiver, maintained by Andreas Müller (not related to me, AFAIK), but I think it's for a long-gone/unsupported version of GNU Radio, so compiling that on a modern system would be a challenge already. However, looking at the source code and flow graphs might be pretty enlightening already! Mr. Michael Höhn implemented he receiver as his Master's thesis, and he came up with this structure (p. 67):

In my head, this diagram pretty neatly translates to a GNU Radio flow graph – that is, a graph built from GNU Radio signal processing blocks. GR is a system to exchange data between these signal processing modules and to implement such yourself. As such, it lends itself very well to implementations of receivers that understand their operation as flow of signal and data through multiple processing stages. In fact, there's a very popular graphical flow graph design tool (the GNU Radio companion) which represents this abstract idea of a signal flow graph as an actual image of connected rectangular blocks. I think, just to get an idea of what SDR implementations might look like, the GNU Radio Guided Tutorials might be a good start.

Furthermore, if you're interested in audio broadcasting, DRM (digital radio mondiale) might be of interest to you, too – and with gr-drm there's a working GNU Radio transmitter implementation for that.

I really do not want to stand between your enthusiasm and diving heads-on into software defined radio and digital signal processing – but you need to be aware that of all engineering disciplines, these are among the more mathematical, algebraically influenced. So your initial "I want to understand the whole machinery of DSP" might be a little too ambitious – I don't think there's anyone on here claiming that they'd do that – it's a huge field of expertise, and everyone can just know the basics and be an expert on their own specialization. There's audio processing, image processing, software defined radio, biosignal, radar, sonar, … experts on here, and some that I look up to as being pretty universally educated in DSP, but none of these is expert in everything. That is a good message for you – you don't have to be an expert in everything to be good in what you do!

So, personally, I might recommend starting slowly, using a Book on Signal Theory, which nowadays always introduces to the basics of DSP, too. Based on that, go for digital communications. You don't have to fully understand every bit of e.g. the information theory used in channel coders and source coders (source coders==compression algorithms), but it's necessary that you understand their role, so having a background in signal theory is not really optional if you want to understand how to recover your data from a very noisy observation.

Have a lot of fun!


"Noise Spectral Density" a.k.a "Power Spectral Density" in datasheets In the case of statistics, there are two terms: "Variance" and "Standard Deviation" to - for in the latter case - have one to note the case in which the units are without squares.

Am I correct to understand, that there is no such differentiation in signal processing, and in both cases - taken root or not - it's always termed "Power Spectral Density"?

I am asking as on page 3 of Analog Devices AD8307 there's the unit of $\frac{nV}{\sqrt{Hz}}$ for the "Input Noise Spectral Density".

Is there a specific reason as to why manufacturers are taking the square-root?

Power spectral density is given in power per $Hz$. "Input Noise Spectral Density" is given in units of volts per $\sqrt{Hz}$.  These are related:

$$
PSD = \frac{W}{Hz}
$$

$$
W = \frac{V^{2}}{R}
$$

$$
Input\ noise\ spectral\ density\ = \frac{V}{\sqrt{Hz}}
$$

The $\sqrt{Hz}$ is from the conversion power to volts.  For example if you want the noise ($nV$) in $1kHz$ $BW$ you must multiply $\frac{nV}{\sqrt{Hz}}\sqrt{1000}$


Typical shape detection algorithms Here is my target, given a set of consecutive discrete 1-d data points within a window. I want to know whether there are good algorithms to detect certain typical shapes of the signal in the window. Including:

A positive step
A negative step
An upward line
A downward line
A hill shape
A bowl shape

The size of the given window could be of any length. And the data is from financial markets so it could be highly random. I wonder if there are good methods to detect whether the above typical shapes exist in the window, which could of course be given by a probability for each shape. I don't need very concrete descriptions of each algorithm, just need to get a direction for my further research.

Do you know Calculus?

All those shapes correspond to bounded values over an interval of the function or its first and second derivatives.

There are many ways to approximate derivatives in discrete sequences.

You can design a metric for each of your conditions based on the derivatives, then evaluate those metrics on a sliding basis.

1) A step occurs when the first derivative is relatively close to zero, then has a spike, then returns to zero.

2) Same as 1, opposite sign.

3) An upward line exists when the first derivative is bound near a positive value and the second derivative is bound near zero.

4) Same as 3, opposite sign

5) A hill (approximated by a parabola) occurs when the second derivative is bound near a negative value.

6) Same as 5, opposite sign


What is the relation between the carrier frequency and the attenuation of electromagnetic waves? For electromagnetic waves, what is the relation between the the carrier frequency, attenuation (absorption loss) and propagation distance?
I know that modulating high frequency carriers are required to reduce antenna size, and to transmit several signals simultaneously.

My question is : does the attenuation increase with these high frequecies when transmitted through air? If so, how is this problem solved?

60GHz is associated with absorption by Oxygen.

Rain and Snow reduce range.

Radio waves can penetrate fresh water.  The Great Lakes was a place where submarine crews trained during WW2 and could send and receive underwater. Salt water is another situation entirely. Frequencies need to be very low.  There is often a good natural wave guide just above the ocean surface.

Foliage can attenuate.

The situation in the optical and IR regions are complicated. There is a program called MODTRAN that models absorption in detail.

For 60GHz, attenuation can actually be a virtue, not a problem. You can pack more receivers and transmitters closer together, at a given complexity.


how do you compute the channel gain from path loss index in wireless communication? how do you compute the channel gain from path loss index in wireless communication?

You will have to properly define path loss index. As soon as you have done that, you'll have a formula that links path loss index and attenuation/gain.


Is DSP a branch of Computational Science? The way I understand it, computational science covers a wide range of fields, all of which use computers and mathematical models to solve a problem involving one of the sciences (physical, life etc.).

But is DSP, which does use maths to manipulate a signal, a form of computational science?

Nope. Its not . Although , computational sciences covers very wide range of subjects but to be able to implement DSP , one must have a thorough understanding of signal fundamentals . That is not a part of computational sciences .


FIR filter to delay signal About filter design, mostly we want to control frequency amplitude response and the linearity of phase response. I wonder whether there are FIR filters specialized to delay signal?

For instance, let's say I want to process a multi-frequencies signal, delaying 1200Hz for offset = $\frac{\pi}{8}$ ($f_s$ = 44100) and leaving 500Hz unchanged.

Is such a FIR filter possible? How can I do it?

It is easily possible to design FIR filters with a complex frequency response specification, i.e., with a desired frequency response

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

where both the desired magnitude $|D(\omega)|$ as well as the desired phase $\phi(\omega)$ can be prescribed. A least squares approximation of $(1)$ by an FIR filter just involves the solution of a system of linear equations. I've written a Matlab/Octave function solving that problem: lslevin.m.

Note that it is important to choose the specification $(1)$ in such a way that it can be realistically achieved by an FIR filter of a given length $N$. E.g., the average delay implied by the desired phase $\phi(\omega)$ should be close to $(N-1)/2$, which is the delay of a linear phase FIR filter of length $N$.


Zero-pad to chirp signal I have generated 16 chirps in experimental study:

and the signal was picked up by some electronic sensor, and then I have calculated the spectrum of the received signal:

Following that I have calculated the cross spectrum density between the sound emitter and receiver and did some post processing, i.e.:

x=voltage1(:,2);
y1=vertgeo1(:,2);

Fs=1000;
window=hamming(512);
noverlap=512/2
nfft=512;Fs=1000;

[Pxy1,f]=cpsd(x,y1,window,noverlap,nfft,Fs);

From this code the number of the points that shows the frequency resolution, essentially the number of the point that f contains is $1\times 256$.

I would like to know how I can increase this to $1025$ points.

No you don't need to zero-pad each chirp. MATLAB's cspd function divide the signal into segments of length nfft. Would you want is to zero-pad each segment. Luckily for you, the cspd function can do it for you. Take a look here:

  If $\tt nfft$ is greater than the signal length, the data is zero-padded. If $\tt nfft$ is less than the signal length, the segment is wrapped so that the length is equal to $\tt nfft$.

So what you should do is increase the value of nfft from $512$ to $2048$.


Period, amplitude and phase through FFT in MATLAB I need to compute period, amplitude and phase of the signal from spectrum analyser. It looks like a sine wave. Since i am new to MATLAB and also have limited idea about Fourier transforms, I am looking for a simplified explanation of how i might use FFT to compute period, phase and amplitude through MATLAB programming. The information i have regarding signal are: sweep time: 0.3183s sweep: 2001 save loop: 100. Thank you so much in anticipation.

Suppose your signal is x.

L = length(x);                   %%% Length of your signal
nfft = L;                        %%% number of points in your fft
res = fft(x,nfft)/nfft;          %%% taking fft of your signal
f = fs/2*linspace(0,1,nfft/2+1); %%% defining correct frequency axis
res = res(1:nfft/2+1);           % half spectrum becos fft result is mirror
figure, plot(f,abs(res));        % here you can amplitude in freq domain
figure, plot(f,angle(res));      % here you can see phase in freq domain.

Hopefully, this will be helpful to you.


Phase of the sine signal by finding FFT of it in matalb using angle I am using sine signal with  phase=0 an frequency=100KHz.
Computed FFT of it.
Using angle function and multiplied the result with 180/pi;
But getting 90.36 at 100KHz frequency.
Can anybody explain why I am getting 90 degree instead of 0 degree?

The basis function for the Fourier Transform is the complex exponential $e^{j\omega t}$
Since
$$sin(\omega t) = \frac{1}{2j}(e^{j\omega t} - e^{-j\omega t}) $$
the Fourier Transform of a sine wave a phase of -90 degrees at the positive frequency and a phase of +90 degrees at the negative frequency.
For a cosine, the phase would be 0 for both the positive and the negative frequency.


How to find the actual size of an existing object from the pixels of an image? I am quite new to vision and OpenCV, so forgive me if this is a stupid question but I have got really confused.

My aim is to detect an object in an image and estimate its actual size. Assume for now I only want length and width, not depth.
Lets say I can detect the object, find its size(length and width) in pixels and I have both the intrinsic and extrinsic parameters of the camera.

The documentation of camera calibration as I understand says that the intrinsic and extrinsic camera parameters can be used to transform from camera coordinate system to world coordinate system. So this should mean converting from pixel coordinates to real coordinates right? And so I should be able to use pixel size and these parameters to find the real size?

But, say the object is photographed at different depth (distance from camera), then its size would come out different using above method. So.....what does it mean transforming from camera to world coordinates?

Any kind of explanation/links/help would be highly appreciated.

I really am no expert in cameras (nor in OpenCV in which I'm pretty sure everything can be handled automatically et efficiently without doing calculations yourself). But without knowing the distance between the lens and your object, that will be difficult.

Actually, knowing the focal length of your camera and the distance from the nodal point to the image plane, you can deduce the distance from the nodal point to the object. Then just use Thales' theorem and you can do exactly what you want.

The problem of this method is that it relies on how well the focus is done on the camera. And because of the varying depth of field of your camera, the magnitude of the error will vary too. AFAIK, this method should work pretty well for short distances. For longer ones, that can be very difficult as the focal point is near infinity and the depth of field is really large.

As an alternative, you can put an object of known size next to the object which size you are measuring. Thales' theorem again and you get your data.


Causal unstable system turn into stable anticausal? I would appreciate it very much if someone would be able to provide some clarity, help or comment on this problem.
I have been reading several papers on time series identification such as

https://www.jstor.org/stable/24305996
https://www.jstor.org/stable/24307363
https://www.jstor.org/stable/24308842

and system identification

https://doi.org/10.1016/S1474-6670(17)56646-X
https://doi.org/10.1109/78.258141
https://doi.org/10.1109/9.827371

among several others.

The question is, given an unstable causal system, is it possible to turn it into an anticausal but stable system in order to simulate it? (assuming that I have the complete input sequence).

I think it might have to do with the time reversal property of the Z-Transform, i.e.

$\mathcal{Z} \left\{ x[n] \right \} = X(z)$ with it's ROC given by $r_2<|z|<r_1$, then

$\mathcal{Z} \left\{ x[-n] \right \} = X(z^{-1})$ with $\frac{1}{r_1}<|z|<\frac{1}{r_2}$

however i can't figure it out.

In the Forssell and Hjalmarsson paper says that backward filtering is common in off-line signal processing applications, like the zero-phase filtering (filtfilt in Matlab) and the implementation of non-causal Wiener filters. The relevant (at least for my question) equations are the (31), (32) and (33).

To illustrate the problem, please consider the following system given by

$y[n] = \frac{B(z)}{A(z)} u[n] = \frac{(z-0.1)}{(z-0.5)(z-5)} u[n]$;

assume that the complete input sequence $\left\{ u[n] \right\}^{N}_{n=0}$ is known (available).

Then, it is possible to separate the system as follows

$y[n] = \frac{B(z)}{A_s(z)} \frac{1}{A_a(z)} u[n],$

where $A_s(z)$ and $A_a(z)$ represents the stable and anti-stable monic parts of $A(z)$ respectively. Thus we can obtain the causal part of $y_{\text{C}}[n]$ by computing

$y_{\text{C}}[n] = \frac{B(z)}{A_s(z)} u[n] = \frac{(z-0.1)}{(z-0.5)} u[n]$

followed by

$y[n] = \frac{1}{A_a(z)} y_{\text{C}}[n] = \frac{1}{(z-5)} y_{\text{C}}[n]$

So here is my problem. The approach given in Forssell and Hjalmarsson is the following

in order to obtain $y[n]$ backward-filter the reversed sequence given by $\left\{ y_{\text{C}}^R[n] \right\}_{n=0}^{N}$ through the filter $\frac{1}{A_a(z)}$ as follows

$y^R[n] = \frac{1}{A_a(z)} y_{\text{C}}^R[n]$

the desired signal $y[n]$ is then finally obtained by reversing the sequence $\left\{ y^R[n] \right\}_{n=0}^{N}$

I apologize in advance for the length of the question and the orthographic mistakes.
Any comments, suggestions and help will be greatly appreciated. Thanks

An unstable system can be stablized by all-pass decomposition for exampe, but simply constructing the anti-causal impulse respone extension will not provide the same output; so I think it's not possible, but I'm not rigorous at this point. Following just shows why.

Assume that your unstable and causal IIR filter has the impulse response $h_+[n]$ and by evaluating it into $ n < 0$ domain you construct a stable but anti-causal impulse response denoted by $h_{-}[n]$.

Your true output given by the unstable filtering is
$$ y[n] = h_{+}[n] \star x[n]$$

Now you construct alternate output from a reversed input and the anti-causal filter as
$$ y_{-}[n] = h_{-}[n] \star x[-n] $$

And by reversing $y_{-}[n]$ , do you get $y[n]$ ?

No as seen:

$$y_{-}[-n] = h_{-}[-n] \star x[n] \neq h_+[n] \star x[n] $$

that's to say unless $h_{-}[-n] = h_+[n]$ , they wont be the same. And indeed the reversed anti-causal impuse response will not be equal to the causal one, since one is stable and the other is not.


Autocorrelation function of a triangular wave I am interested to find the analytical expression for the autocorrelation function of a signal that comprises triangular pulses. I have followed the derivation in "Statistical Theory of Communication" by Lee (1960, p.16)
However, if I replicate this in MATLAB, then I do not obtain the autocorrelation function as derived here. The general shape is the same, but it looks as follows:

with the autocorrelation function dipping below zero. In this case, the triangular pulses are spaced such that each triangle is half of the period, with the remainder of each period being zero.
Why is this different from Lee's Derivation, and how can I go about adapting his technique to find the analytical form that matches the empirical result?
The code I used in MATLAB was as follows:
function [x, acf1, acf2] = ACFTriangle(leftbase,rightbase,height,gap,repeat)

x = [zeros(1,gap) (1:leftbase)*(height/leftbase) (rightbase-1:-1:1)*(height/rightbase) 0];

x = repmat(x,[1,repeat]);

t = 0:length(x)-1;

acf = xcov(x,'unbiased');

n = length(acf);

s = (n-1)/2 + 1;

acf1 = acf(s:n)/var(x);

acf2 = autocorr(x,length(x)-1);

x = [t; x]';

end

The reason for my question is that I have a very large number of time series in which there may be a hidden triangular pulse wave. If I have an analytical expression for the autocorrelation sequence, then I can identify their characteristics from the shape of the autocorrelation function.
New edit
Notwithstanding the excellent answer below, how would one go about modifying Lee's mathematical derivation to allow for removal of the means? I can attempt this doing the following:
\begin{equation}
\varphi_{11} = \frac{1}{T_1}\int_0^{b-\tau} \left(\frac{E_m}{b}t - m\right)\left(\frac{E_m}{b}\left(t+\tau\right) - m\right)\,dt
\end{equation}
where $m$ is the mean of the series.  However, following this through gives:
\begin{equation}
\varphi_{11} = \frac{E_m^2}{6b^2T_1}\left(\tau^3-3b^2\tau + 2b^3\right) + m^2\left(b-\tau\right) - \frac{mE_m}{bT_1}\left(b-\tau\right)^2 - \frac{mE_m}{bT_1}\tau\left(b-\tau\right)
\end{equation}
which clearly gives $\varphi_{11}=0$ for $\tau=b$.  Where is the error in this working?

I couldn't quite follow your code (you calculate two different versions of the ACF?), but I believe the problem with the plot being shifted toward zero is that xcov calculates the cross covariance not correlation. Remember that cross covariance is subtracting the mean, so there should be a shift! You should be calling xcorr instead.
The following code snippet reproduces the Figure 9 from your screenshot.
tau = 1;
b = 100;
Em = sqrt(6);
T1 = 2*b;
signal = zeros(T1, 1);
step = Em/(b+1);
signal(tau:tau+b) = linspace(0, Em, b+1);
signal = repmat(signal, 1000, 1);
[xcorrOut, lag] = xcorr(signal, 'unbiased');

figure
hold on
grid on
plot(lag, xcorrOut)
xlim([-5*b, 5*b])

Edit
In response to the new edit in the question, the OP wants a way to convert the cross correlation result to a cross covariance result. The cross covariance can be obtained from the cross correlation by subtracting the product of the mean and conjugate mean, https://en.wikipedia.org/wiki/Cross-covariance:
$$ \varphi_{11}(\tau) - |\mu_{f_1}|^2$$
where $\mu_{f_1}$ is the mean of $f_1(t)$, as shown in Figure 7 of the screen capture. This additional code snippet builds upon the above to show how adjusting the correlation can get you the same result as computing the covariance.
mu = mean(signal);
xcovFromXcorr = xcorrOut - mu*conj(mu);
[xcovOut, lag] = xcov(signal, 'unbiased');

figure
hold on
grid on
plot(lag, xcovFromXcorr, 'LineWidth', 2)
plot(lag, xcovOut, '--', 'LineWidth', 2)
xlim([-5*b, 5*b])


Nonnegative or positive band-limited interpolation Given samples of an everywhere non-negative or positive-valued continuous-time signal band-limited to half the sampling frequency, is there some practically applicable way to interpolate it so that there is no danger of producing negative values? The approach should be scalable to arbitrary accuracy so basic linear interpolation won't do. By arbitrary accuracy I mean an arbitrary small difference between the band-limited continuous-time signal and the interpolation.

Perfect sinc interpolation would work because negative contributions would get cancelled by positive contributions. Ideal sinc interpolation is not possible in practice so other filters are used. This could give negative values because usually the filter impulse response has negative values. For example cubic Hermite spline piece-wise polynomial interpolation gives this impulse response:

Figure 1. Cubic Hermite spline impulse response has negative values.

Lagrange polynomials used as splines would also have negative values in the impulse response. B-splines have everywhere non-negative impulse responses, but their pass-band frequency responses deviate increasingly from the brick-wall frequency response of sinc interpolation with increasing B-spline order. The B-spline impulse response tends to a Gaussian function as the order is increased. A Gaussian function in time domain is a Gaussian function in frequency domain, so not a perfect brick-wall filter.

The obvious way to fix the problem would be to 1) take the absolute value of the interpolation output or 2) set all negative output values to zero, but I'm wondering if there is some better (inherently non-negative) method, because both 1) and 2) would create discontinuities in the derivative of the signal. This also creates spectral images.

I don't know how bad the problem stated in this question really is, because if the band-limited signal is non-negative, this also rules out some contradictory discrete signals like $[\dots, 0, 0, 0, 1, 0, 0, 0, \dots].$ The corresponding non-negative continuous signal would then have to be sinc, which is partially negative, which is the contradiction.

Shay Maymon has investigated "bandlimited square-roots" in his 2011 thesis Sampling and Quantization for Optimal Reconstruction. The idea there is to find (not easy) another band-limited signal (complex and at half the bandwidth) the squared magnitude of which equals the wanted band-limited signal. Interpolating that other signal and multiplying the result by its complex conjugate is never negative. That reminds me of this question about interpolation of the magnitude of discrete Fourier transform (DFT).

If we require the impulse response of the interpolation filter to be non-negative then we can obviously guarantee that the interpolated signal is also non-negative. However, it can be shown that the frequency response of a low pass filter with non-negative impulse response is always far from ideal. More specifically, it can be shown that if $\epsilon$ is the maximum negative pass band deviation of the squared magnitude $|H(\omega)|^2$ of the filter's transfer function from its value at $\omega=0$

$$|H(\omega)|^2\ge |H(0)|^2-\epsilon,\quad |\omega|<\omega_c \tag{1}$$

where $\omega_c$ is the cut-off frequency, then the following inequality holds:

$$|H(\omega)|^2\ge |H(0)|^2-4\epsilon,\quad |\omega|<2\omega_c\tag{2}$$

This inequality imposes a severe restriction on the minimum stop band attenuation if the ripple size $\epsilon$ is chosen sufficiently small (as should be the case for a good low pass filter).

As an example, if $\epsilon$ is chosen such that the pass band ripple is $1$ dB, the attenuation at $\omega=2\omega_c$ cannot be better than $7.5$ dB.

Of course, if $\epsilon\ge |H(0)|^2/4$ then $(2)$ represents no restriction, but in that case the maximum pass band deviation of the magnitude of the transfer function is greater than usually required for a good low pass filter.

Since the inequality $(2)$ is independent of the filter order, an interpolation filter with arbitrary accuracy cannot be realized by a filter with non-negative impulse response.

The proof of $(2)$ can be found in Attenuation Limits for Filters with Monotonic Step Response by A. Papoulis, IRE Transactions on Circuit Theory ( Volume: 9, Issue: 1, March 1962).


How do I get texture feature images from GLCM texture analysis? In this paper under section 2. "Proposed Methodology", the author creates a GLCM from an image and then extracts texture features resulting in a new image for each feature.  I've seen other authors do this in other papers as well.  The texture features used in the image are the common Haralick features detailed on page 619 here.

How did the author extract a texture feature image from those features?

The author has a unique image for contrast, an image for entropy, an image for correlation, etc.  When I look at the Haralick feature equations the result is always a number not another image.
How is a unique image created for each texture feature?

EDIT:

For GLCM definition as suggested in comments, a GLCM stands for the Gray-level co-occurence matrix.  It's a common technique used for texture analysis.  It is defined over an image to be the distribution of co-occurring pixel values at a given offset.

Looking at these references, this is what I would understand:

The Haralick descriptors only create a single value for each GLCM (as you also thought).
The Rampun paper creates 32 images, describing 32 different features for each input image.
The crucial point is this sentence: On the other hand, we used
a small window size of 5 x 5 throughout the process
This sentence, in combination with this one (taken from here, p.25, Sec. 2.6): The GLCM Texture Image is the result of moving the GLCM window across the entire image or sub-image.

So, my understanding is:

Create moving window of size 5x5
move this window over your input image
for each window position (x,y), calculate a feature (e.g. contrast), and write the value for this feature into the "contrast" image at position (x,y).

This way, you extract "local" features, i.e. the GLCM is only done over local area of size 5x5. This also makes more sence, since the Rampun paper wants to segment images, i.e. it needs local features instead of global ones.


Preparations to implement FFT on a wav file (code Java, not Library) I am making preparations to implement FFT on a wav File.wav Only Java.
I'm reluctant to use matlab

I have a File wav stereo, 16 bits.

I need to calculate the FFT, but I have confusing.

My File has a Frequency of 44100, the samples in bytes are 10 minuts ->

(10 minuts)*(60 seconds)*(44100 sample rate)*(2 bytes because are 16 bits)
= 52920000 >~ 52MB per channel.

As you can see it's a big file.

1. Question.
Can I to split my file in order to make smaller calculus?
If yes, My file can be variable, How I can calculate a moderate lenght working in function of frequency (greater frequency will need greater length FFT).

2. Question.
According to FFT the size need to have a lenght = mantissa = Pow(2, Integer).
Can I to fill with zeros to reach Pow(2, Integer), Example I have a file with a lenght of 65000 samples, then I will need a FFT with order of 16 -> mantissa = 65536, Can i to fill with zeros the others 536 samples?

Thank you.

Simply taking the Fourier transform of a ten minute signal is unlikely to provide you with any useful information - I suspect what is actually of interest to you is the spectral information at intervals throughout the signal, e.g. 600 frequency plots at 1 second intervals, which if you so desire can be combined into a spectographic image for analysis. This can be achieved by performing multiple windowed Short Time Fourier Transforms (STFT) and storing the results.

Can I to split my file in order to make smaller [calculations]?

Yes, it is up to you what granularity of frequency information you wish to obtain. Small FFT lengths will afford you a better time resolution, longer FFT lengths will afford you a better frequency resolution. For purely visualisation purposes an FFT length of 1024 is generally sufficient.

Can I to fill with zeros to reach Pow(2, Integer) [for my FFT length]?

Yes, if you do not have sufficient samples to reach a power of two length for your FFT then you should zero pad the input. This will introduce the minimal amount of spectral error into your frequency domain results.


I have a graph where I need to find a particular point based on sampling I currently have a graph which looks something like this:

Now I have discrete data that represents the Number of Clusters(the x values), giving us the Average Cluster Distance (the y values). There is only a single elbow point, after which the values just descend to 0. Finally, I do know the maximum amount that K can take (the x domain). Based on this I have three questions:

Is there a way I can find the value of the elbow point through differentials (first or second order)?
In which case, what is the smallest sampling rate that will allow me to get a graph with the desired characteristics? I don't care too much about the accuracy as of right now (though I would like to know some kind of a tool that lets me measure the error). If I cannot use differentials, I am still interested in understanding how I could recreate the graph with the lowest sampling rate, and thus find the elbow point with the least amount of processing (For E.g. Nyquist Theorem, I would have 1/2 the amount of elements and still recreate the signal).
Finally, would incorporating random selection into the selection of the signal help me in any way?

Based on the graph you posted and the elbow defined, all you have to do is to take the first difference of the signal and compare it to a threshold value, assuming that the signal is monotonically decreasing. Such as this:

if( (y[k] - y[k-1]) < Th )
          elbow = k-1

if the signal is not monotonically decreasing, you should look for other measures as well.

Sorry but I cannot tell the required sampling rate for this graph as you stated it comes from an already discrete data...?


Is spectral density conserved after aliasing? Will the integral of the PSD of a down-sampled signal still equal the variance in all cases

No. The aliased component will interfere with the non-aliased components and the interference can constructive or destructive.
Trivial example:
$$x[n] = \sin\left(\frac\pi2n\right)$$
If you down sample this to $y[n] = x[2n]$, you get all zeros.


The essential bandwidth of a rectangular pulse The problem is to estimate the essential bandwidth of a rectangular pulse

\begin{equation}
g(t) = \Pi(t/T),
\end{equation}

Where the essential bandwidth must contain at least $90\%$ of the pulse energy.

Here's what I'm trying to solve :

$$0.9 =\frac1\pi\int\limits_0^{2BT\pi} \mathrm{sinc}^2 \left(\frac x2\right) \,dx$$

I'm not sure how to get MATLAB to solve this numerically for the upper bound value.

first the energy for the pulse is $E_g$

\begin{equation}
E_g = \int^{+\infty}_{-\infty} g^2(t)dt\, = \int^{+T/2}_{-T/2} dt\, = T
\end{equation}

and

\begin{equation}
\ \Pi(t/T) \Longleftrightarrow T\mathrm{sinc}(\pi fT)
\end{equation}

and for the Energy Spectral Density of the pulse

\begin{equation}
\ \Psi(f) = \mid G(f) \mid ^2 = T^2\mathrm{sinc}^2(\pi fT)
\end{equation}

The energy $E_B$ within the band from $0$ to $B$ hertz is

\begin{equation}
\ E_B = \int^{+B}_{-B} T^2 \mathrm{sinc}^2(\pi ft)df\,
\end{equation}

setting

\begin{equation}
\ 2 \pi fT = x
\end{equation}

\begin{equation}
\ df = dx/((2 \pi fT)
\end{equation}
plugging into the integral
\begin{equation}
\  E_B = (T/ \pi)  \int_{0}^{2 \pi BT} \mathrm{sinc}^2(x/2)dx\,
\end{equation}
\begin{equation}
\ E_g = T
\end{equation}

\begin{equation}
\ 0.9 = (E_b/E_g) =  1/ \pi \int^{2 \pi BT}_{0} \mathrm{sinc}^2(x/2)dx\,
\end{equation}

I wanted to work through this example fully, but I'm not sure how to get this numerically solved, solving for the bandwidth B from the upper bound of the integral using MATLAB.

First note that $$\mathcal{F}\{\mathrm{rect}(\frac{t}{T})\}=T\mathrm{sinc}(Tf)=\frac{\sin(\pi Tf)}{\pi f}$$
Then you have correctly calculated the pulse energy as $$E_g=\int_{-\infty}^{+\infty}g^2(t)dt=T$$
Now we want to find $B$ such that
$$\int_{-B}^{B}T^2\mathrm{sinc}^2(Tf)df=0.9T$$
or since the function is even,
$$\int_{0}^{B}T\mathrm{sinc}^2(Tf)df=0.45$$ which by assuming $$Tf=x\Rightarrow df=dx/T$$ is
$$\int_{0}^{TB}\mathrm{sinc}^2(x)dx=0.45$$
Using the following code we get the essential bandwidth $$\bbox[0.5em,#efe,border:0.1em groove navy]{\
B\approx\frac{0.8455}{T}\
}$$

J = 0.45; %target value
acceptable_error = 1e-4;
dx = 1e-5;
xmax = 1; %initial guess
max_iterations = 100;
iteration = 0;
while iteration < max_iterations;
    x = 0:dx:xmax;
    y=sinc(x).^2;
    I = trapz(x,y);
    if abs(I-J) <= acceptable_error
        break;
    end
    if I < J
        xmax = 3/2*xmax;
    else
        xmax = xmax/2;
    end
    iteration = iteration+1;
end
fprintf('absolute error = %f\n',abs(I-J));
fprintf('solution = %f\n',xmax)


All-pass filters to remedy group delay I would like to sum a filterbank made of bandpass butterworth filters to reconstruct an audio signal. Given group delays for each channel shown in the following plot . . .

How would I go about designing complementaty allpass filters to flatten this out somewhat? It doesn't have to be completely flat, just a bit flatter than it is. Sampling the impulse responses to make linear phase filters is not an option at the moment.

The Matlab code used to generate this is shown below. Note: I have multiplied the group delays by a factor of 2 as the signal is passed through the bank twice.

sr = 22050;

%Make octave wide filters
cfs = [125 250 500 1000 2000 4000];
bw = 1; %octaves

nBins = 2^12;
gd = zeros(nBins,numel(cfs));

for nn = 1:numel(cfs)
    cf = cfs(nn);
    loEdge = cf*(2^-bw/2);
    hiEdge = cf*(2^+bw/2);
    [b, a] = butter(2, [loEdge hiEdge]/(sr/2), 'bandpass');
    [gd(:,nn),f] = grpdelay(b,a,nBins,sr);
end

figure; semilogx(f,2*gd); xlim([10 8000])
legend(num2str(cfs')); xlabel ('F [Hz]'); ylabel('Group delay [samples]')

Maybe you'd be better off using Bessel-Thompson filters. Maybe one could find an all-pass that approximately corrects a Butterworth filter, but it looks rather hard.


MATLAB: $\tt freqz$ vs $\tt bode$ I started studying DSP recently. In my MATLAB assignments, I am often asked to plot the frequency response of a digital filter using freqz. I know there's another function for plotting a frequency response, bode.

In my last assignment, I was asked to design a notch filter for $60\textrm{ Hz}$ with a sampling frequency of $200\textrm{ Hz}$. After doing so, I plotted the frequency response using freqz, which gave me the expected result - the notch was on $60\textrm{ Hz}$:

However, when I plot the frequency response using bode, the frequency is off significantly:

To my understanding, both functions measure the frequency response. Why is one giving me the correct result then, when the other is not?

In fact both plots are the same.

Your freqz plot uses a linear discrete frequency, $\omega$ axis scaled between $0$ and $1$ with units of $\pi$ rad per sample. For 200 Hz sampling frequency, a 60 Hz notch would be shown at $w_c = 0.6 \pi$ rads, which is put at 0.6.

The bode plot, however, uses a logarithmic scale for the continuous-time radian frequency $\Omega$. Given that you have a $f_c=60$ Hz of notch frequency, its radian value is $\Omega_c = 2\pi \times  60 \approx 377$ radians per second.

So both plots are showing the notch at the correct place.


find the gradient direction at a point in a shape I am trying to implement the generalized Hough transform in matlab. The algorithm requires the gradient direction at each point in the shape. How can I measure $ \phi $ as shown in the figure below?

Assuming you are using discrete inputs and that s(0..n-1) (or s(1...n) in Matlab) is an array of x,y points along the boundary of a 2d shape ...

Gradient(i) = (s(i).x, s(i).y) - (s(i-1).x, s(i-1).y)
Amplitude(i) = sqrt(Gradient(i).x^2 + Gradient(i).y^2)
Angle(i) = atan(Gradient(i).y / Gradient(i).x)

To make you calculations relative to (Xc,Yc) substract (Xc,Yc) from each point in s prior to calculations.

If your calculations are continuous, replace the above s(i)-s(i-1) with limit(s)/limit(sigma) and use the symbolic toolkit.


How to name every row of 100*2500 matrix in MATLAB? currently i am working with a matrix. The dimension of the matrix is 100*2500. I need to assign each row the matrix to a new variable. How i can do this?

You can assign each row to a vector as follows:

row1 = matrix(1,:);

That will store the whole first row of matrix in row1.

To do this "automatically", you can store all the vectors containing each row in a cell array, for example:

for k = 1:size(matrix,1) %For each row
    rows{k} = matrix(k,:); %Store the k-th row in the k-th element of the cell array
end


What is the cause of slow roll off of windowed filters Exactly why using the of windowing function, such as Hamming or Blackman, causes a slower transition to stop band? Everywhere it is mentioned just as an obvious fact that this happens but I haven't found a text that mentions what is the cause.

An example for what I mean:

An FIR low-pass filter designed with cut off frequency 1 Hz without a window function reaches stop band at little over this frequency (1.04 Hz or so).
When window function is applied, the filter reaches the stop band at over 1.2 Hz

Also, is there a way to compensate for this? For example, if I really want the cut off frequency to be 1 Hz and use window function, should I design the filter to have cut off say 0.9 Hz?

A (non-full-width-rectangular) window attenuates data at its edges, thus narrowing the width (decreasing the 1st moment) of the bulk hump of the input.  Decreasing width in the time domain increases uncertainty in the frequency domain, which is represented by a filter response which has a shallower transition band stretching farther from the center frequency.

One way to roughly compensate is to start with a wider data window (longer FIR), which when "narrowed" by a non-rectangular window function will constitute roughly the same bulk "width" as before the widening.


Should the input of a Kalman filter always be a signal and its derivative? I always see the Kalman filter used with such input data. For example, the inputs are commonly a position and the correspondent velocity:

$$
(x, \dfrac{dx}{dt})
$$

In my case, I only have 2D positions and angles at each sample time:

$$
P_i(x_i, y_i) \qquad \text{and} \qquad (\alpha_1, \alpha_2, \alpha_3)
$$

Should I compute velocities for each point and for each angle to be able to fit the Kalman framework?

A state variable and its derivative are often included as inputs to a Kalman filter, but this is not required. The essence of the Kalman framework is that the system in question has some internal state that you are trying to estimate. You estimate those state variables based on your measurements of that system's observables over time. In many cases, you can't directly measure the state that you're interested in estimating, but if you know a relationship between your measurements and the internal state variables, you can use the Kalman framework for your problem.

There is a good example of this on the Wikipedia page. In that example, 1-dimensional linear motion of an object is considered. The object's state variables consist of its position versus time and its velocity on the one-dimensional line of movement. The example assumes that the only observable is the object's position versus time; its velocity is not observed directly. Therefore, the filter structure "infers" the velocity estimate based on the position measurements and the known relationship between velocity and position (i.e. $\dot{x_k} \approx \frac{(x_k - x_{k-1})}{\Delta t}$ if acceleration is assumed to be slowly-varying).