Analyze and reproduce sonic screwdriver sound Do you know about Doctor Who and its screwdriver?

Well, I'm trying to understand how to replicate this sound but the spectral analysis is too way complicated, just see its spectrogram.

I was trying to figure out what kind of sine waves or partials it is made of, but probably I have no good experience in this.

I'm writing it in python but for know I have nothing than a signal with two sine waves just to make tests.

How can I understand how this sound is made?
There is a frequency modulation component which I would not consider right now, because it is probably a steady state sound which was modulated later.

import matplotlib.pyplot as plt
import numpy as np
import sys
from scipy.signal import get_window
import sounddevice as sd
import soundfile as sf

timeLength = 1.0                       # seconds
fs = 44100
t = np.arange(0, timeLength, 1.0/fs)
ww = get_window("hann", t.size)
A = .01
phi = np.pi # phase, radiants.

T = 1.0/fs
df = 20
t_circ = 0.1
t_mid = 0.2
f_mid = 3700
f_circ = df*(1.0 - ((t_circ - t_mid)/t_mid)**2.0)**0.5
phase1 = np.cumsum((f_mid + f_circ)*T*np.ones(t.size))
phase2 = np.cumsum((f_mid - f_circ)*T*np.ones(t.size))
circ = np.sin(2.0*np.pi*phase1) + np.sin(2.0*np.pi*phase2)

ss = []
#brown_noise = np.random.wald(1, 0.002, len(t))
for i in range(10):
    f = f_mid + 20*i
    s = A * np.cos(2 * np.pi * f * t + circ)
    ss.append(s)
ss = sum(ss)
ss = ss/max(ss)

plt.plot(ss)
plt.show()
sf.write('out.wav', ss, fs)

An approximation

Actually, the signal is not that complicated in my opinion. However, I am not a sounds guy... Just playing around with it a bit.

Interestingly, I have a very similar signal that I have produced (see attached sound file and spectrograms. Spectrogram above is my signal (somewhat tuned to your signal's parameters) and below is the sonic screw-driver (of which Doctor by the way? ;) )

I cannot post the code here, since as I said the parameters are a bit tuned to your signal. What I basically did is to create a signal that is...

A sum of 10 signals. Center frequency f_mid = 3.7 kHz + 20Hz*N, N = 0, 1, ..., 9, just to have some kind of a wobbly jitter in there
The center frequency is modulated by a "circle" (df: maximum change of the center frequency, t_circ: duration of a circle, t_mid: starting point of the circle in time) with df = 20 Hz, t_circ = 0.1s. The phase of the signal is generated from the constant frequency plus the dynamic modulation and a very simple numerical integration (phase multiplied by an increasing dt-vector and cumsuming over it). To get to a signal's frequency, you would differentiate its phase wrt to time - so this is the inverse process.

T = 1.0/44100.0
f_circ = df*(1.0 - ((t_circ - t_mid)/t_mid)**2.0)**0.5
phase1 = np.cumsum((f_mid + f_circ)Tnp.ones(t_circ.shape))
phase2 = np.cumsum((f_mid - f_circ)Tnp.ones(t_circ.shape))
circ = np.sin(2.0*pi*phase1) + np.sin(2.0*pi*phase2)

Added some somewhat "Brownian" Noise, however it is much below the noise in your signal.

Things to explore

The "duration" of the "circles" seems so vary slightly. Probably sinusoidal variation is also fine, so the period could be jittered.
Apparently, there is additionally a very broad noise-band. Maybe try some uniform noise with a time-varying windowing?
Subtract a bunch of short-term, constant frequency signals from your signal. There is a nearly some kind of a grid visible in the spectrogram.

Some code samples

I am sorry, in my previous description, I apparently have missed some details. Please find here a code snipped that should be able to produce a basic signal that exhibits the basic "circular" frequency change:

import matplotlib
matplotlib.use('MacOSX')
import matplotlib.pyplot as pp
import numpy as np
import scipy.io.wavfile
import scipy.signal

# Export function for convenience
def export_wave2(t, signal_l, signal_r, name = "test_sp.wav"):
    sample_rate = 44100
    data = np.stack((signal_l, signal_r))
    data = np.transpose(data)
    scipy.io.wavfile.write(name, sample_rate, data.astype("float32"))

# Sampling frequency
fs = 44100
T = 1.0/fs

# Circle parameters
t_circ_max = 0.2 # duration of one "circle"
t_circ     = np.arange(0, t_circ_max, 1.0/fs)
t_mid      = t_circ.max()/2.0
f_mid      = 3700 # middle frequency of the sound signal
df         = 400     # amplitude of the frequency modulation, ie frequency will vary between f_mid - fd and f_mid + fd

# Create ONE circle
f_circ = df*(1.0 - ((t_circ - t_mid)/t_mid)**2.0)**0.5

phase1 = np.cumsum((f_mid + f_circ)*T*np.ones(t_circ.size))
phase2 = np.cumsum((f_mid - f_circ)*T*np.ones(t_circ.size))
circ = np.sin(2.0*np.pi*phase1) + np.sin(2.0*np.pi*phase2)

# Repeate the circle N times, ie the signal will be N*t_circ_max long
N = 10
circs = np.hstack([circ for k in range(0,N)])
circs = circs / np.abs(circs).max()
circs_time = np.arange(0, N*t_circ_max, 1/fs)

pp.plot(circs_time, circs)
pp.show()

export_wave2(circs_time, circs, circs, name = "test_sp.wav")

The frequency behaviour is of course not visible directly in the plots, but when a spectrogram is calculated from circs. I have loaded the test_sp.wav into Audacity for this.

Sound example


FFT - radix 2 - decimation in time -frequency range
Hello,

I would like to ask how to decrease (make it narrow) frequency range for calculations in FFT radix 2 - decimation in time - algorithm?

Let me explain what I exactly mean:

When you compute DFT in regular manner (I mean not FFT) you make frequency bin loop, and for each frequency bin you need next loop to use each possible sample you have. Then you can measure whole frequency range - in example for audio it would be something like from 0 Hz to 44100 Hz, of course if you have buffer size the same as sample rate 44100. But very often we don’t need to calculate whole frequency range, in example for audio very often 20 Hz to 20 kHz is perfectly enough. So in regular computation of DFT, you can just change frequency bin loop range, and then you use all possible samples but measure only those chosen frequences.

But when you compute DFT by FFT algorithm, you first divide your buffer size by half of a buffer size, so in my example you compute 22050 DFTs with two point size of each. And then you conquer results in log2(buffer size) steps. From the the first step to the one before last there is no loop that correspond to frequency range or just I can’t see it. Only the last step correspond to frequency range.

So my question is:

can I shrink the frequency range only in that last step, or can I do that in some way for whole FFT algorithm? If yes, could you give me some hint, how to do that?

If you studied the materials I recommended you on comments to your previous question. You would notice that FFT algorithm takes advantage of certain symmetries to make computations faster. If you don't want to work with a whole frequency range simply select the appropriate sampling frequency and concentrate all your FFT resolution points there. If you can't choose your sampling frequency take the samples and then decimate  them. Alternatively you can also use the Goertzel algorithm to work only with the frequencies you select.


How significant is the time lag between inertial IMU data in a inertial analisys? I'm working on a MCU system that use an IMU (ICM-20948) which has two different chip inside: one provides accelerometer and gyroscope data at $1125$ Hz (configurable at $125$ Hz), the other one magnetometer data at fixed $100$ Hz (not configurable).
Now I'm wondering how much the non-synchronicity of the data could be a problem in a motion processing algorithms.

It depends on your algorithm, and on whether the timing is known.
Assuming that the timing of each measurement is known, you can incorporate measurements at unsynchronized moments by executing a Kalman update step for the measurements you have, at the moment the measurements happened -- just use a measurement matrix ($H$ in most Kalman literature) that matches the measurement you have, and do your update.
This does mean that you need to use the appropriate state update equations, which makes for the filter design being a royal pain in the behind -- but it's more a matter of dotting a bunch of 'i's and crossing a bunch of 't's (and not crossing 'i's or dotting 't's) than it is a matter of the basic math getting worse.


What is the meaning of this notation? I have seen this notation in the book Schaums-Outline a couple of times now and I can't quite understand it. I'd very much appreciate if someone explained what it means. I am thinking that it means Transform, as in the input goes in the system and transforms? Thank you very much!

It just means "the transformation that turns $x$ into $y$."

You might also see $\mathbf{T}^{-1}$ which means the inverse: turning $y$ into $x$.


What is Super Resolution in the Context of Windowing? I have been trying to understand what super-resolution is, in the context of DSP/DIP. What criteria is being maximized/minimized, and why?

Most of my online searches yield super-res techniques from the optical physics stand point, however I know anecdotally that they have been used in image processing, and in some doppler-radar processing as well.

One particular instance of Super-Res I remember seeing was in the design of a windowing function, where the main lobe was the same width as that of a boxcar, while simultaneously having very low and fast degrading side lobes, similar to a hamming window. In this sense, the 'super-res' came from the fact that frequencies very close to one another were able to be resolved due to the low main lobe width, and a very high dynamic range. Is this what super-res really amounts to?

Is there an explanation that, (perhaps), likens it to a form of de-noising, or sparsity pursuit? Or is it a different animal altogether?

The term super-resolution is used very loosely nowadays. But I think the following problem was the original idea, or at least, the one that made the term famous.
Suppose you have a scene and several observations of that scene, i.e. several frames of a video with slight camera motion between them. Super-resolution algorithms combine those several observation into one highly resolved image. This is possible because each observation might contain unique information not present in the others.
If you consider a certain frame as the key frame, this process can be interpreted as a MAP estimation since a priori (or extra) information is present in others frames.
Additionally, we can use natural image statistics as extra information to regularize the solution. Total variation regularization is a broadly used technique.
I suggest the paper Sung Cheol Park - Super Resolution Image Reconstruction: A Technical Overview.


Image watermarking in wavelet domain (data hiding in images) I am working and researching on Image watermarking/steganography and I'm stuck at a point and want to know which technique is better in image processing for data hiding.

Bit Plane Complexity Segmentation OR CDMA Spread spectrum.

Can some experts give their opinion on both techniques and which one is better and why? like advantages and disadvantages or comparison of both

Thanks a lot.

So, I just want to hide some text in a color or grayscale image using wavelet trasform.

I could use either bpcs or cdma spread spectrum technique. I want to know which one should go for and why.

This is for learning purpose and no particular application.

Thanks. Please let me know if u need any specific info.

For image watermarking, I suggest you not to use wavelet Transform.
Instead, you can use an efficient technique in terms of robustness and transparency, based on the Quantization Index Modulation (QIM). You can also use a variant of it which is called "Spread Transform Dither Modulation".
Good luck


image correlation without fft I want to stitch 2 images by crosscorrelation measure, but I don't want to use fft.

I tryed this but something seems to be wrong.any suggestions?

void test(IplImage* src,IplImage* temp)
{
double max_peak=0;
                Point pt;
                double per_x=0.5;
                double per_y=0.5;
                for(int cy=-h*per_y;cy<h*per_y;++cy)
                {
                        for(int cx=-w*per_x;cx<w*per_x;++cx)
                        {
                                uchar* s_p= (uchar*) src->imageData;
                                uchar* t_p= (uchar*) temp->imageData;
                                __int64 suma=0;
                                __int64 sumb=0;
                                for(int y=0;y<h;++y)
                                {
                                        for(int x=0;x<w;++x)
                                        {
                                                suma+= s_p[x];
                                                if((y+cy)<h&&(x+cx)<w&&
                                                        (y+cy)>0&&(x+cx)>0)
                                                        sumb+= t_p[x+cx];
                                        }
                                        s_p+=w;
                                        t_p+=w;
                                }
                                int meana = (int)(suma / (w * h)),
                                        meanb = (int)(sumb / (w * h));

                                s_p= (uchar*) src->imageData;
                                t_p= (uchar*) temp->imageData;

                                __int64 corr = 0,
                                        denoma = 0,
                                        denomb = 0;

                                for(int y=0;y<h;++y)
                                {
                                        for(int x=0;x<w;++x)
                                        {
                                                int da = (int)s_p[x] - meana;
                                                int db= 0;
                                                if((y+cy)<h&&(x+cx)<w&&
                                                        (y+cy)>0&&(x+cx)>0)
                                                        db = (int)t_p[x+cx] - meanb;
                                                else
                                                        db= -meanb;

                                                corr += da * db;
                                                denoma += da * da;
                                                denomb += db * db;
                                        }
                                        s_p+=w;
                                        t_p+=w;
                                }

                                double peak= corr / sqrt( (double)denoma * denomb);
                                if (peak>max_peak)
                                {
                                        pt.x= cx;
                                        pt.y= cy;
                                        max_peak=peak;
                                }
                        }
                }
}

At the start of the inner loop, you assign

uchar* t_p= (uchar*) temp->imageData;

later on, you use it as

t_p[x+cx];

Aren't you missing the cy here?


IIR Cascaded Filter Response, Unexpected Results I am attempting to understand how to manipulate portions of a cascaded series of biquad notch filters in a control system I am working with.

As part of determining the composite response of all the cascaded filters, I am starting with looking at the individual biquad filters. I have taken the existing control system coefficients that I have for the filters and have saved the various biquad filters in the Filter design tool. When I use the MATLAB/Simulink filter design tool, I obtain an unexpected response as well as an unstable root/pole plot. I have excerpted an image below of one of the filters from the control system.

The system sampling frequency is 10 khz. The Form is Direct Form II. The a0 coefficient is normalized to 1.

My specific questions are:

Is this response still consistent with a notch filter?
If this is still applicable as a notch filter, why would the filter start with a negative magnitude and then return to a positive magnitude? For context, I am used to the analog realm and a notch filter starting at 0 db and ending at 0db.
If I needed to attenuate a specific frequency as viewed from power spectra of the feedback in the control system, would it be appropriate to modify a particular notch filter of my cascaded series to focus in on that frequency or since they are cascaded, would I have to modify all terms slightly?

Thanks for any insight.

Looks like a simple sign error in the denominator. If you use a = [1 -1.866685 .91058], you get something that looks notch-like at around 380 Hz.

  The Form is Direct Form II

This is numerically the worst. Shouldn't matter here, but it's good practice to either stick with Direct Form I or Transposed Form II.


Impulse response of a transfer function Lets say we have a transfer function from an LTI-system that goes as follows, and we want the impulse response for it:
$$\frac{10\cdot e^{\left(-\mathrm{i}\right)\cdot \pi\cdot f\cdot 8}}{5+\mathrm{i}\cdot 2\cdot \pi\cdot f}$$
The way I would approach solving this, would be to use the following fourier transform formulas:
$$\frac{1}{a+i2\pi*f}=exp(-at)*u(t)=\frac{10}{5+i2\pi*f}=10*exp(-5t)*u(t)$$
$$G(f)*exp(-i*2*\pi*f*t_0)=g(t-t_0)$$
$$exp(-i*\pi*f*8)=exp(-i*\pi*2*f*4) =>t_0=4$$
$$exp(-i*2*\pi*f*4)=g(t-4)$$
Which would end up with the following result:
$$\frac{10\cdot e^{\left(-\mathrm{i}\right)\cdot \pi\cdot f\cdot 8}}{5+\mathrm{i}\cdot 2\cdot \pi\cdot f}=10*exp(-5t)*u(t)*g(t-4)$$
However, this seems to be incorrect. Did I do an error with the fourier transforms, or have I misunderstood some of the formulas?

You're almost there; you just need to connect a few dots.  Let your $\frac{10}{5 + i 2 \pi f} = G(f)$.  Then $g(t) = 10 u(t) e^{-5 t}$.

Now we get into that exponent part.  Your $h(t) = g(t - t_0)$ is correct, but you're applying it incorrectly.  You need to apply it to the part that I've labeled as $g(t)$: $h(t) = g(t - t_0) = 10 u(t - t_0) e^{-5 (t - t_0)} = 10 u(t - 4) e^{-5 (t - 4)}$.

So, basically, it's the output you'd expect from the "pure" transfer function, just delayed by 4 time units (presumably seconds).

Does that help?


How to find the group delay of IIR filters in Python/SciPy I sent a square wave signal (red) through a IIR filter (butterworth in this case), to make it more "realistic". Then there's a slight group delay in my output signal (blue). Is there a way to find out what the group delay is? I only need the total group delay, not spectrum of group delay. I need the number so I can shift my decision point accordingly, so I am still sampling at the middle of the eye diagram to achieve optimum result.

I only need the total group delay, not spectrum of group delay.

Group delay is a spectrum, so this doesn't make sense.  The group delay is the derivative of the phase response of the filter, so in Python it can be calculated as

from scipy import signal
from numpy import pi, diff, unwrap, angle

w, h = signal.freqs(b, a)
group_delay = -diff(unwrap(angle(h))) / diff(w)

The units of phase are radians, and the units of frequency are radians/sec, so the units of group delay are in seconds.  Examples at 1 and 2

The group delay of a Bessel filter is maximally flat, so when you say "the" group delay, you probably mean the maximum passband group delay.  Like in this filter:

The passband has about 0.032 seconds of group delay.  Your Butterworth filter does not attempt to have a flat group delay in the passband, however, so it's up to you to decide what "the" group delay is:

Do you want the group delay at 0 Hz, or the maximum group delay at the cutoff frequency?  Or maybe at a specific frequency like the fundamental clock frequency of your digital signal?


Relationship between input and output sequence in Hartley transformation As you know that Discrete Hartley transformation is related to the discrete Fourier transformation, $i.e$, assuming we have a vector $X = [x_0,x_1,\ldots,x_N]$, its Hartley transformation is equal to $H(X) = Real(FFT(X)) - Imag(FFT(X))$. where $H$ denotes the Hartley transformation.

I wonder, if I want the output of the Hartley transformation equals to a vector of length $N$ whose all elements are $1$, it means $H(X) = [1,1,\ldots,1] $, what supposed to be the input $X$? I means I need to formulate the relationship between the inputs and outputs.

Wikipedia's entry for the discrete Hartley transform shows states that the $\mathsf{DHT}$ is, up to a scaling, its own inverse. If $x$ is a vector with $N$ entries and $y$ is its discrete Hartley transform,
\begin{equation}
y = \mathsf{DHT}x,
\end{equation}
then
\begin{equation}
x = \frac{1}{N}\mathsf{DHT}y.
\end{equation}

If $x$ is a vector with $N$ entries such that
\begin{equation}
\mathsf{DHT}x = \underbrace{(1,1,\ldots,1)^{\mathsf{T}}}_{\textrm{$N$ entries}},
\end{equation}
then we recover $x$ with
\begin{equation}
x = \frac{1}{N}\mathsf{DHT}\left(\begin{array}{c}1\\1\\\vdots\\1\end{array}\right).
\end{equation}
This means that $x$ is
\begin{equation}
\begin{split}
x &=~ \frac{1}{N}\left(\mathsf{Re}\left[\mathsf{DFT}\left(\begin{array}{c}1\\1\\\vdots\\1\end{array}\right)\right] - \mathsf{Im}\left[\mathsf{DFT}\left(\begin{array}{c}1\\1\\\vdots\\1\end{array}\right)\right]\right),
\end{split}
\end{equation}
where $\mathsf{DFT}$ is the discrete Fourier transform, which we usually compute with a FFT algorithm. The $\ell^{\textrm{th}}$ entry of the $\mathsf{DFT}$ of the all-1 vector is
\begin{equation}
\begin{split}
\sum_{n=0}^{N-1}1\times e^{-2\pi j \ell n/N}
&=~ 1 + e^{-2\pi j \ell/N} + \left(e^{-2\pi j \ell/N}\right)^2 + \cdots + \left(e^{-2\pi j \ell/N}\right)^{N-1}\\
&=~ \left\{\begin{array}{rl}N & \textrm{if}~\ell=0,\\0&\textrm{if}~\ell\neq 0.\end{array}\right.\\
&=~ N\delta_{\ell,0},
\end{split}
\end{equation}
where $\delta_{p,q}$ is the Kronecker delta. One way to show this is to note that if $\ell = 0$, then each exponent is $0$, so each term in the sum is $1$. On the other hand, if $\ell\neq 0$, then $\exp(-2\pi j\ell/N) \neq 1$,and
\begin{equation}
\begin{split}
1 + e^{-2\pi j \ell/N} + \left(e^{-2\pi j \ell/N}\right)^2 + \cdots + \left(e^{-2\pi j \ell/N}\right)^{N-1}
&=~ \frac{1 - \left(e^{-2\pi j \ell/N}\right)^N}{1 - e^{-2\pi j \ell/N}}\\
&=~ \frac{1 - e^{-2N\pi j \ell/N}}{1 - e^{-2\pi j \ell/N}}\\
&=~ \frac{1 - 1}{1 - e^{-2\pi j \ell/N}} ~=~ 0.
\end{split}
\end{equation}

That shows that the $\mathsf{DFT}$ of the all-1 vector has no imaginary part, and its real part is $(N,0,0,\ldots,0)^{\mathsf{T}}$. Hence
\begin{equation}
x ~=~ \frac{1}{N}\left(\begin{array}{c}
N\\0\\0\\\vdots\\0
\end{array}\right)
~=~ \left(\begin{array}{c}
1\\0\\0\\\vdots\\0
\end{array}\right).
\end{equation}


Spatial sampling frequency equation Upon reading this microphone array tutorial, I'm stuck in the following piece:

The spatial sampling frequency along the x-axis is given by
$f_{x_s} = \frac{\sin \theta \cos \phi}{\lambda}$

where $\theta$ and $\phi$ are the elevation and azimuth angles of the incoming signal, and $\lambda$ is the wavelength.
I cannot grasp the meaning of the equation. How come the sampling frequency depends on the angles? Also, shouldn't it depend on the inter-element distance?

The spatial sampling frequency of a linear array with equidistant sensor (antenna) spacing on the $x$-axis is simply given by

$$f_x=\frac{1}{d}\tag{1}$$

where $d$ is the distance between the sensors (cf. Eq. $(35)$ in the document you linked to). What is meant by the equation in your question is the spatial frequency of the incoming wave along the $x$-axis, according to the $x$-component of the direction vector $\mathbf{\alpha}$ (cf. Eq. $(13)$ in the document).

To avoid spatial aliasing you need to guarantee

$$\frac{1}{d}>2\mathbf{\alpha}_x=2\frac{\sin\theta\cos\phi}{\lambda}\tag{2}$$

for all possible angles $\theta$ and $\phi$. The maximum value of the numerator of the fraction on the right-hand side of $(2)$ is $1$, and the minimum value of the denominator is the minimum possible wavelength $\lambda_{\text{min}}$:

$$\frac{\sin\theta\cos\phi}{\lambda}\le \frac{1}{\lambda_{\text{min}}}\tag{3}$$

Combining $(2)$ and $(3)$ gives

$$d<\frac{\lambda_{\text{min}}}{2}\tag{4}$$

which corresponds to Eq. $(38)$ in the document.


how to eliminate frame boundary discontinuity? I'm applying at each 2048 samples the pitch synchronous algorithm to pitch scale some audio files, everything seems to occur well, but I've noticed some glitches, artifacts, clicks between every boundary frame!

I'm testing off-line in large áudio files(I'm processing frame-by-frame), what I tried to do was keep the old, current and the next frame to pitch synchronous computation (then I pitch scale 2048*3 samples), but does not seem to be working as expected.

What are the known ways to eliminate frame boundary discontinuity?

Windowing is the first line of defense.  you may need the windowing to be complementary, which is that all window functions (at their various frame delays) add to 1.  the Hann window is an example.  the Hamming window (which might have better sidelobe behavior than the Hann) is sorta complementary, but it adds to 1.08 (just a simple scaling issue).

if your window use is only for analysis (and not for synthesis) you might not need it to be complementary, so you could use a Kaiser window or a Gaussian window, both which have some useful mathematical properties in analysis.

the second line of defense is to use your 2048-sample as just buffers and space your frames (and the analysis windows that come with), by the detected period and adaptively advance the frame locations by a period each frame process.  that's a little weirder because you need to have, a priori, some knowledge of the period which what your pitch synchronicity is all about.

the most weirdest line of defense, is to advance each frame by a well-chosen "random" value between $L$ and $2L$ (where $2L$ is longer than the longest anticipated period) to scramble up whatever effect of the boundary discontinuity.

as we might say on comp.dsp, "you'll have to pay for more information".


CTFT to DTFT why can't you always just substitute $\Omega = \omega/T_S$ This is something I've always wondered about in DSP class, but just accept as a fact because I never really fully understand why this is the case:

Given CTFT:

$$X_s(j\Omega) = 6000 \pi \sum \limits_{k=-\infty}^{\infty} \Bigg[ \delta(\Omega - 1200\pi k - 4000 \pi)
 ~+~  \delta(\Omega - 1200\pi k + 4000 \pi)\Bigg]$$

Then converting to DTFT using formula:

$$X(e^{j\omega})=\frac{1}{T_s} X_s\bigg(\frac{\omega}{T_s}\bigg)$$

But, Why is answer this:

$$X(e^{j\omega}) = 6000 \pi \sum \limits_{k=-\infty}^{\infty} \Bigg[ \delta(\omega + T_s(-1200\pi k - 4000 \pi))
 ~+~  \delta(\omega +T_s(- 1200\pi k + 4000 \pi))\Bigg]$$

instead of this:

$$X(e^{j\omega}) = 6000 \pi \sum \limits_{k=-\infty}^{\infty} \Bigg[ \delta(\frac{\omega}{T_s} - 1200\pi k - 4000 \pi)
 ~+~  \delta(\frac{\omega}{T_s} - 1200\pi k + 4000 \pi)\Bigg]$$

seems like this formula: $$X(e^{j\omega}) = \frac{1}{T_s} X_s\bigg( \frac{\omega}{T_s}\bigg)$$ should really be:  $$X(e^{j\omega}) = \frac{1}{T_s} X_s\bigg( \omega = \Omega T_s\bigg)$$

It is due to the property of unit impulse function

\begin{equation}
\delta(\alpha t)  = \frac{ 1}{|\alpha|} \delta(t)
\end{equation}

(Usually $\omega $ is used for angular frequency and $\Omega $ is used for angle. Interpret my answer accordingly. )

suppose
\begin{equation}
X_a(\omega) = \delta(\omega -x_0 )
\end{equation}
and if we substitute '$\omega$' with $\frac{\Omega}{T_s}$ then
\begin{equation}
X_a(\frac{\Omega}{T_s}) = \delta(\frac{\Omega}{T_s}  -x_0 )
   = \delta(\frac{\Omega - x_0 T_s}{T_s} )
   =  \delta(\frac{ 1}{T_s} (\Omega - x_0 T_s) )
   =  {T_s} *\delta(\Omega - x_0 T_s)
\end{equation}
where '$\alpha$'  is $\frac{ 1}{T_s}$
and
$$
\frac{ 1}{|\alpha|} = T_s
$$

and when you do DTFT calculation on

\begin{equation}
X(\Omega)  = \frac{ 1}{T_s} X_a(\frac{\Omega}{T_s})
   =\frac{ 1}{T_s} * ( T_s * \delta(\Omega - x_0 T_s) )
   = \delta(\Omega - x_0 T_s)
\end{equation}

in the final answer $T_s$ and $\frac{ 1}{T_s}$ cancels

and you will get the expression based on the independent variable '$\Omega - x_0 T_s$'

So there is nothing wrong in your actual answer.


Compensate for phase shift of Hilbert transform I'm using Octave and trying to use the "instfreq" function from the time frequency toolbox which requires an analytic signal input. To get this analytic signal I'm using the Octave function "hilbert" on my real valued signal thus:

analytic_sig = hilbert( my_real_valued_signal ) ;
[ instf, t ] = instfreq( analytic_signal ) ;

where instf is the instantaneous frequency I'm trying to obtain. However, the results are not quite what I expected. I suspect that the 90 degree phase shift induced by the hilbert transform is where the problem lies, so my question is how can I alter the output of the Octave "hilbert" function to compensate/adjust for this phase shift?

Edit in response to comments

With a simple test case I'm creating a signal x with instantaneous frequency ifl

octave:1> [x ifl] = fmlin(10,0.05,0.35)
x =

  -0.97815 - 0.20791i
  -0.80902 - 0.58779i
  -0.30902 - 0.95106i
   0.50000 - 0.86603i
   1.00000 + 0.00000i
   0.30902 + 0.95106i
  -0.91355 + 0.40674i
  -0.30902 - 0.95106i
   1.00000 - 0.00000i
  -0.50000 + 0.86603i

ifl =

   0.050000
   0.083333
   0.116667
   0.150000
   0.183333
   0.216667
   0.250000
   0.283333
   0.316667
   0.350000

but in my envisioned real life application I will only have data for the real component

octave:2> real_sig = real(x)
real_sig =

 -0.97815
 -0.80902
 -0.30902
  0.50000
  1.00000
  0.30902
 -0.91355
 -0.30902
  1.00000
 -0.50000

using the hilbert function I create an analytic signal from my available real data

octave:3> anal_sig = hilbert( real_sig )
anal_sig =

  -0.97815 + 0.07265i
  -0.80902 - 0.41187i
  -0.30902 - 0.92331i
   0.50000 - 0.81514i
   1.00000 + 0.04490i
   0.30902 + 0.98765i
  -0.91355 + 0.52573i
  -0.30902 - 0.89042i
   1.00000 + 0.28002i
  -0.50000 + 1.12978i

and when this is put in to the instfreq function I get the measured instantaneous frequency instf

octave:4> instf = instfreq( anal_sig )
instf =

   0.10520
   0.13131
   0.15427
   0.18209
   0.20487
   0.24755
   0.31328
   0.30974

The instfreq function returns values in the range [2:end-1] of its input, so padding to make an easy comparison

octave:5> [ ifl [ 0 ; instf ; 0 ] ]
ans =

   0.05000   0.00000
   0.08333   0.10520
   0.11667   0.13131
   0.15000   0.15427
   0.18333   0.18209
   0.21667   0.20487
   0.25000   0.24755
   0.28333   0.31328
   0.31667   0.30974
   0.35000   0.00000

it can be seen that the measured frequency is different from the true, known frequency. Thinking of a phasor diagram, I was thinking that perhaps some trigonometric manipulation of the real and imaginary components of the hilbert function output might correct this.

1) A Hilbert transform has a very very long impulse response (above some given noise floor), so you need a ton more data to manufacture an analytic signal, otherwise you won't have enough to span the width of the Hilbert impulse response filter without serious edge truncation effects.

2) Instantaneous frequency estimates from this type of artificially approximated analytic signal can be very noisy, so often has to be low pass filtered to significantly below half your sample rate.


Help me understanding the super equalizer source code Please forgive my bad English.

Hello, I'm trying to understand Shibatch's super equalizer Winamp plugin source,
but I'm not good at DSP and mathematics.

It seems that super equalizer plugin uses overlap-add method to evaluate FIR filter.
The code that calculates the FIR filter coefficients in super equalizer is like this:
(Just the simplified pseudocode... not the exact code.)

for (int i = 0; i != FILTER_LENGTH; ++i)
{
    int t = i - FILTER_LENGTH / 2;

    result[i] =
        band_gain[0] * LPF(t, band_freq[0]) +
        band_gain[1] * (LPF(t, band_freq[1]) - LPF(a, band_freq[0])) +
        band_gain[2] * (LPF(t, band_freq[2]) - LPF(a, band_freq[1])) +
        band_gain[3] * (LPF(t, band_freq[3]) - LPF(a, band_freq[2])) +
        band_gain[4] * (IMPULSE(t) - LPF(t, band_freq[3]));

    result[i] *= WINDOW(t);
}

LPF is lowpass filter function. LPF(t, CutoffFrequency)

IMPULSE(t) is 1 when t is 0. Otherwise IMPULSE(t) is 0.

WINDOW(t) is Kaiser window function.

band_freq array contains equalizer frequency bands.

band_gain array contains the weight of each frequency bands.

I assumed there are just four frequency bands.

band_gain[0] = weight of 0Hz ~ band_freq[0]Hz
band_gain[1] = weight of band_freq[0]Hz ~ band_freq[1]Hz
band_gain[2] = weight of band_freq[1]Hz ~ band_freq[2]Hz
band_gain[3] = weight of band_freq[2]Hz ~ band_freq[3]Hz
band_gain[4] = weight of band_freq[3]Hz ~ SAMPLING_RATE / 2

I want to know this part:

result[i] =
    band_gain[0] * LPF(t, band_freq[0]) +
    band_gain[1] * (LPF(t, band_freq[1]) - LPF(a, band_freq[0])) +
    band_gain[2] * (LPF(t, band_freq[2]) - LPF(a, band_freq[1])) +
    band_gain[3] * (LPF(t, band_freq[3]) - LPF(a, band_freq[2])) +
    band_gain[4] * (IMPULSE(t) - LPF(t, band_freq[3]));

What formula does the code using?
Or... what should I study to understand the code?

Any answer is welcome.
Keyword suggestion is welcome.

Thanks.

result[i] =
band_gain[0] * LPF(t, band_freq[0]) +
band_gain[1] * (LPF(t, band_freq[1]) - LPF(a, band_freq[0])) +
band_gain[2] * (LPF(t, band_freq[2]) - LPF(a, band_freq[1])) +
band_gain[3] * (LPF(t, band_freq[3]) - LPF(a, band_freq[2])) +
band_gain[4] * (IMPULSE(t) - LPF(t, band_freq[3]));

is actually very simple to explain with pictures. First we assume f0 < f1 < f2 < f3.... lets look at one section

band_gain[0] * LPF(t, band_freq[0]) means create a low pass filter with cutoff frequency band_freq[0] it then says to multiply the amplitude of the filter by band_gain[0] something like this

now lets move onto the more complicated lines

band_gain[1] * (LPF(t, band_freq[1]) - LPF(a, band_freq[0])) means make a LPF with Cutoff band_freq[1] then subtract from it a LPF with cutoff band_freq[0] if you look at the image below thats the pink section. Then you multiply the result by `band_gain1'

you can find the other bands in a similar fashion . Now lets add them together, something like this, just 'superposition'

which makes our final filter look roughly like this

I hope those pictures make sense


Image zooming with Fourier transform I want to go from this image

into this one:

So basically I need to scale the white square. The authors of the paper claim that this can be done in four steps:

zero-padding in real space (image is treated as 1D vector).
discrete Fourier transform.
cropping the central part of the spectrum in frequency domain.
inverse Fourier transform.

I don't see how it is being done in practice and if this works. Can you help out?

What I believe is being done is interpolation in the frequency domain to increase the number of pixels while smoothing the transition from black to white.

Consider a single row for example, which is a rectangular function in time. The Fourier Transform would be a Sinc function, which would dominate the outer part of the spectrum in the frequency domain (go toward zero in the middle, which represents $F_s/2$ where $F_s$ is the sampling rate, and then rise back as a mirror of the first half up to $F_s$. Adding zeros in time simply interpolates more samples of this same function in frequency. Then truncating the center of the spectrum reduces the higher frequency range of the Fourier Transform, which effectively makes the pulse wider (lower frequency).

I discuss zero padding for FFT's in more detail here: What happens when N increases in N-point DFT


Why upsample before modulation? We're doing a project in which we're sending an OFDM frame, modulated to some appropriate carrier, over a channel. Before modulating, we are instructed to upsample and low-pass it. Is there a good explanation this should be done? We have full control over sample rates at both Tx and Rx.

I guess it has something to do with an interpolated signal resembling a smooth signal more, which makes it somehow more appropriate for multiplying by a modulating function. But this is very imprecise and hand-wavy.

Can anyone explain why one would want to upsample, given that there are no differences in sample rate between transmitter/receiver?

The reason is aliasing. Because the modulation is digital, on the discrete signal, the sampling frequency should be high enough that the channel's frequency band does not extend beyond the Nyquist frequency = sampling frequency / 2. If the sampling frequency is too low, the frequencies beyond the Nyquist frequency will alias to other frequencies, meaning that you would be transmitting at the wrong frequency. By lowpass filtering before modulation, you are removing spectral images of the baseband. If those are not removed the signal has frequencies outside the baseband and similarly the modulated signal won't be contained in the channel's band.


Signal response in Frequency Hz from angular frequency I have a filter whose respose I have in angular frequency.

A0 = [1.0000000,0.3281250,0.0000000];
B0 = [0.5859375,-0.0625000,0.0000000];
[z,p,k]=tf2zp(B0,A0);
figure;
zplane(z,p);
w = 0:(pi/500):pi;
H=freqz(B0,A0,w)
figure;
plot(w,20*log10(abs(H)/max(abs(H))))

1) How can I get the corresponding response in Hertz? Is it:

plot(w/(2*pi),20*log10(abs(H)/max(abs(H))))

2) Also, how do I get the response of the above filter, if I apply it to a sinusoidal signal at 20MHz, with a sampling freq of Fs=300MHz

have a look to this example: http://www.mathworks.com/matlabcentral/fileexchange/55822-low-pass-filtering

here you have a signal that can be represented/plot in both frequencies scales:

omegas_axis=[0:2*pi/L:2*pi-2*pi/L];%axis input signal
frec_axis=[0:fs/L:fs-fs/L];%axis input signal

the "omegas_axis" is the angluar frequency.

Hope this helps,

BTS.


Cross Power Spectral Density of Unevenly Sampled Data Here's my problem. The input signals $x$ and $y$ will be having the time value aligned with each other. However, the data are not evenly sampled. I would like to calculate CPSD of both signals.

The solution comes to my mind is as follow

$R_{xy}$ = cross correlation of $x$ and $y$ (I'm not sure how to do it with unevenly sampled data)
Use Lomb Periodogram to estimate the CPSD of Rxy

Is there anything I miss? And what is the right way to do it?

Update

Sam Maloney suggested to use interpolation to fill the gap to produce evenly sampled data.
This solution is good.

However, the data I obtained may have some missing data / long gap. Interpolating the data may "predict" the data wrongly and produce undesired end result.

As quoted in Numerical Recipes (Spectral Analysis of Unevenly Sampled
Data)

  However, the experience of practitioners of such interpolation techniques is not reassuring. Generally speaking, such techniques perform poorly. Long gaps in the data, for example,
  often produce a spurious bulge of power at low frequencies (wavelengths comparable to gaps)

This is the reason why I choose Lomb Periodogram over FFT in Step 2 without interpolating the data. However, I'm unable to figure out how to solve the Step 1. Is there any way to calculate the cross correlation of 2 signals without interpolating them?

One way would be to interpolate the signals to produce evenly-spaced samples and then calculate the CPSD as normal.


Clarification on the application of window pre-FFT Good morning, I have a simple question about the application of a Hanning window before performing an FFT.  My application is pretty textbook, i'm decoding a 4-FSK signal.

Question: why does the spectral width increase when I apply a Hanning window to a perfect sin?

Many thanks!
A

Generally, any window (rectangular or Von Hann) convolves the input with the transform of the window.

The transform of a rectangular window is a Sinc shaped "ringy" waveform with infinite ripples, but these side lobes are invisible for pure sinusoids that are exactly integer periodic in the FFT width.

Smoother windows have lower ripples, but in exchange have wider main lobes.


Measuring with OpenCV and OpenNI I'm playing with the built in OpenNI access within OpenCV 2.4.0 and I'm trying to measure the distance between two points in the depth map. I've tried this so far:

#include "opencv2/core/core.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"

#include <iostream>

using namespace cv;
using namespace std;

Point startPt(0,0);
Point endPt(0,0);

void onMouse( int event, int x, int y, int flags, void* )
{
    if( event == CV_EVENT_LBUTTONUP) startPt = Point(x,y);
    if( event == CV_EVENT_RBUTTONUP) endPt   = Point(x,y);
}

int main( int argc, char* argv[] ){
    VideoCapture capture;
    capture.open(CV_CAP_OPENNI);
    //capture.set( CV_CAP_PROP_OPENNI_REGISTRATION , 0);
    unsigned t0=clock();
    if( !capture.isOpened() ){
        cout << "Can not open a capture object." << endl;
        return -1;
    }
    unsigned elapsed=clock()-t0;
    cout << "initialized in "<< elapsed <<" s. ready!" << endl;
    namedWindow( "depth", 1 );
    setMouseCallback( "depth", onMouse, 0 );
    for(;;){
        Mat depthMap;
        if( !capture.grab() ){
            cout << "Can not grab images." << endl;
            return -1;
        }else{
            Mat show,real,valid;
            if( capture.retrieve( depthMap, CV_CAP_OPENNI_DEPTH_MAP ) ){
                depthMap.convertTo( show, CV_8UC1, 0.05f);
            }
            capture.retrieve(valid,CV_CAP_OPENNI_VALID_DEPTH_MASK);
            if( capture.retrieve(real, CV_CAP_OPENNI_POINT_CLOUD_MAP)){
                unsigned int sp = valid.at<unsigned char>(startPt.x, startPt.y);
                unsigned int ep = valid.at<unsigned char>(endPt.x, endPt.y);
                if(sp == 255 && ep == 255){
                    Vec3f s = real.at<Vec3f>(startPt.x, startPt.y);
                    Vec3f e = real.at<Vec3f>(endPt.x, endPt.y);
                    float dx = e[0]-s[0];
                    float dy = e[1]-s[1];
                    float dz = e[2]-s[2];
                    float dist = sqrt(dx*dx + dy*dy + dz*dz);
                    putText(show,format("distance: %f m\n",dist),Point(10,10),FONT_HERSHEY_PLAIN,1,Scalar(255));
                }
            }
            circle(show,startPt,3,Scalar(255),3);
            circle(show,endPt,3,Scalar(192),3);
            line(show,startPt,endPt,Scalar(128));
            imshow("depth",show);
        }
        if( waitKey( 30 ) >= 0 )    break;
    }

}

but I have a few issues:

Unless I check for values in the valid depth mask first the program
crashes sometimes (I imagine due to bad data)
As far as I understand the point cloud map returns XYZ data in meters (CV_32FC3), but the numbers I see on screen look wrong.

Am I retrieving the xyz values and doing the distance calculations correctly ? What am I missing out ?

Coordinates should be passed in y,x order, not x,y.
Full updated code:

#include "opencv2/core/core.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"

#include <iostream>

using namespace cv;
using namespace std;

Mat world;
Point startPt(0,0);
Point endPt(0,0);
float dist;

void onMouse( int event, int x, int y, int flags, void* )
{
    if( event == CV_EVENT_LBUTTONUP) startPt = Point(x,y);
    if( event == CV_EVENT_RBUTTONUP) {
        endPt   = Point(x,y);
        Vec3f s = world.at<Vec3f>(startPt.y, startPt.x);
        Vec3f e = world.at<Vec3f>(endPt.y, endPt.x);
        float dx = e[0]-s[0];
        float dy = e[1]-s[1];
        float dz = e[2]-s[2];
        dist = sqrt(dx*dx + dy*dy + dz*dz);
    }
}

int main( /*int argc, char* argv[]*/ ){
    VideoCapture capture;
    capture.open(CV_CAP_OPENNI);
    capture.set( CV_CAP_PROP_OPENNI_REGISTRATION , 0);
    unsigned t0=clock();
    if( !capture.isOpened() ){
        cout << "Can not open a capture object." << endl;
        return -1;
    }
    unsigned elapsed=clock()-t0;
    cout << "initialized in "<< elapsed <<" s. ready!" << endl;
    namedWindow( "depth", 1 );
    setMouseCallback( "depth", onMouse, 0 );
    for(;;){
        if( !capture.grab() ){
            cout << "Can not grab images." << endl;
            return -1;
        }else{
            Mat depthMap,show;
            capture.retrieve(world, CV_CAP_OPENNI_POINT_CLOUD_MAP);
            if( capture.retrieve( depthMap, CV_CAP_OPENNI_DEPTH_MAP ) ) depthMap.convertTo( show, CV_8UC1, 0.05f);
            line(show,startPt,endPt,Scalar(255));
            putText(show,format("distance: %f m",dist),Point(5,15),FONT_HERSHEY_PLAIN,1,Scalar(255));
            imshow("depth",show);
        }
        if( waitKey( 30 ) >= 0 )    break;
    }
}


Pulse detection against various noise floors I'm decoding a pulse amplitude modulation encoded digital signal. The goal is to detect the individual pulse against a background noise floor. The amplitude of the pulses varies depending on the signal strength and quality.

My sampling rate is 2x the bit width such that one bit fits into sample. Thus a single bit is the large value amongst the noise floor (8-bit unsigned values):

6, 1, 18, 100, 9, 4, 7

Detecting the bit is simply looking for the large value relative to its neighbors (100 in the above example). This would be easy if the noise and pulses had known constant thresholds, but in real life they vary. What I actually see is varying noise and pulse amplitude depending on the quality of the signal. Consider these samples (47 and 52 are the bits I'm trying to detect):

2, 8, 12, 47, 19, 13, 10
23, 11, 30, 52, 17, 9, 11

I've been struggling with developing an algorithm that can dynamically determine the threshold and thus pick-out the bits. I'm approaching this with a general computer science background with little to no DSP knowledge.

What are some software techniques to reliability detect a pulse against varying background noise?

I cannot comment on your question, so I am writing an answer here instead.

Be aware that any time you mention 'threshold', you are implicitly making a statement that you are comfortable with those three probabilities:

Certain probability of false alarms, (PFA).
Certain probability of detection, (PD).
Certain probability of misses, (PM).

Put another way, every time you set a threshold, those three numbers will change accordingly. This is inescapable, in any detection problem. Typically one can set the probability of false-alarm that they are willing to 'live with', and go from there.

For your problem, you mention changing noise variance. One way to circumvent this is to Z-Standardize your data vector before attempting anything. This then centers your data, (de-meaning), and forces the variance to unity. This take care of the problem of you having to re-compute thresholds all the time.

With this pre-processing step out of the way, you can now focus on your Test-Statistic, and your Threshold. For example, a test-statistic can simply be each sample. Or, it can be the average of three samples. Or the average of 10 samples. Either way, you will be comparing your test-statistic, against your threshold. Again, you pick a threshold you "can live with", that is to say, one that gives you false-alarms at a rate that your application deems acceptable.

For example, if a doctor has a false alarm rate of 10% in informing you that you have a cold, that might be acceptable for this application, because it is not the end of the world, and the worst thing that can happen is you drink more orange juice than usual. If on the other hand you are designing a fraud detection system for your bank and you false alarm in detecting fraud 10% of the time, you are going to have some very, very angry customers swearing at your representatives all day for locking their credit cards.

Back to your application: It seems that those bits are the result of a demodulated PAM receiver system, so you probably already have a significant amount of averaging gain. This means that a test-statistic of one sample is good enough to go off of. All you need to do now, is set your threshold (of your standardized data), such that you false alarm as some rate you are comfortable with.

(Remember, that you are dealing with standardized data, so you need only set this threshold once, and then forget about it. As a practical note to assist you in determining the threshold, try to plot the histogram of your data without signal, after you have standardized it. This will be your 'noise PDF'. This can assist you in estimating the optimal threshold ala Neyman_Pearson detector.)


How did they color old black and white movies? This question is on dsp.SE as I'm mostly interested in the signal processing part.

There is an Indian movie Mughal-e-Azam which was released in 1960 in black & white which has been reproduced in color in 2004.

How did they color each pixel perfectly?
What technique have they
used to identify the color placement on each pixel?

Look at one of the screenshots from the movie:

I've got an Einstein black-and-white photo which I want to colorize. How is it possible to do so  without knowing what he was wearing back then and what was the actual color of his clothes, background etc.

There is no way of recovering the original color information from a black and white photo, so whether Einstein (resp. Waheeda Rehman) was wearing a pink or green sweater (resp. Dupatta) is up to your imagination.

Historically, this has been done by hand, by painting over the film. The first digital techniques to automate the process consisted in "painting" a few dots of color on each frame, at the center of each uniformly colored region, and using something like a voronoi partition + some blurring to get a color map for each frame (see for example US patent 4606625).

Today, this can be done relatively easily (though manually) with video editing software, by using vector masks to indicate regions of uniform color on a few keyframes, and interpolating between them. Then a color transform is applied to each mask. See it in action here.

Standard image segmentation and region tracking techniques can be used to automate the task of segmentation and the marking of regions on each keyframe - for example by propagating manual annotations to similar/adjacent pixels in space/time, or by detecting uniformly textured regions. Texture and gray level similarity can be used to propagate color cues from a color image to a greyscale image depicting a similar subject - in this case the manual process only consists in finding a template color image - this later task can be itself automated using content-based image retrieval techniques.


Independence of Functions of random Variable Consider I am given two functions of one random variable each for example x=cos(at),y=rect(bt) where a  and b are random variables.And I am given Probability density function for a and b then if I am asked if the two functions are independent or  not so, I want to confirm that before proceeding I will have to convert the pdf of a and b to pdf of x and y or can I directly proceed with pdf of a and b by proceeding I mean checking if joint probability density is equal to product of marginal probability density.

It happens that if $X$ and $Y$ are independent then so will their functions $g(X)$ and $h(Y)$ be; but not $g(X,Y)$ and $h(X,Y)$.


Summing magnitudes of individual frequency bins I am working on signal analyser and have problem with understanding relationship between display resolution and frequency resolution (range).
My frequency range is $ 1-20000\quad [ H z ] $ so to project each frequencies magnitude on display (horizontally in linear scale) I need display horizontal resolution at least 20000.
But of course I don't have :)

So let's say I have 1000 pixels horizontal resolution. So to project whole frequency range I need to express $ \frac { 20000} { 1000} = 20 $ frequency bin magnitudes by one pixel.

And my question is what is the most relevant manner to do that?

Should I just sum all 20 magnitudes? Or better make some avarage value like $ \frac { 1} { 20} \sum _ { i = 1} ^ { 20} mag_{fbin} $

I am pretty sure it's concern a lot of other subjects like: is magnitude print in decibels (log) scale or as range from 0 to 1? Also important thing is perception of human ear (if it's audio). And other things. But I don't know how to deal with that. Could you give any advice? To make the answer easier, let's say I want to create audio signal analyser.

Please help. Thanks in advance.

Having an interface that limits what can be  presented is a common problem. The solution depends a lot on the application and if you want a zoom function.

You can average, take a max, take a median, interpolate, decimate, or whatever seems to work.

Its a cognition problem.


What is leakage in frequency domain adaptive filters? What does it mean by leakage in case of digital filters? My specific question is about the frequency domain adaptive filter function provided in the Matlab DSP toolkit, accessed as adaptfilt.fdaf. It has a parameter called LEAKAGE, but I am not sure what exactly does it represent or how it affects the filter response.
The filter created as h = adaptfilt.fdaf can be used as
[y, e] = filter(h, x, d)
which filters the data in x. Studying the source code of the filter function provided shows how leakage has been used, but the theory on which this filter function is based does not include the leakage component. What does it do and why is it there, if the paper on which this function is based does not have it?

In adaptive filtering, leakage is a stabilization method which may be useful if the covariance matrix is close to singular (i.e. at least one of the eigenvalues is very small), or if there are finite-precision effects in the implementation of the adaptive filter. Leakage changes the update formula such that not only the mean squared error but also the norm of the filter taps is minimized. This prevents unbounded growth of the filter coefficients in cases of numerical ill-conditioning.

For you this simply means that you initially use a leakage factor of 1 (i.e. no leakage) if the FDAF works properly with your input signals. If you encounter coefficient drift (large fluctuation about the optimum solution), you can start by slightly decreasing the leakage factor until the coefficient fluctuation becomes sufficiently small. Note that leakage achieves stabilization at the expense of performance degradation because due to the changed update formula it introduces some bias in the filter taps.


Number of samples in a discrete signals after sampling a continuous periodic signal should be integer? If a continuous periodic signal is sampled with time interval Ts , and T being time period of continuous signal --
For the resultant discrete signal to be periodic following must be true--
the ratio of T/Ts = rational number = number of samples in the discrete signal.

If I take an example of sine wave with T = 300 seconds , Ts = 100 seconds
Number of samples = 3

If I take an example with same sine wave  T = 300 seconds , Ts = 175 seconds
Number of samples = 300 / 175 = 1.714285714285..(repeating)

My question is how come the number of samples in a discrete signals is not an integer , in other words I am not able to understand how can a discrete signal have number of samples in decimals.

There is nothing unordinary about it. Signals can have a period that is not an integer multiple of samples, and that is normal. For example, if you record a 440 Hz sine wave with a sampling rate of 48 kHz, it is not an integer multiple. It takes 11 full sine wave periods to hit an integer amount of 1200 samples.


Number segmentation techniques Assume that the speed limit at a place is captured and is provided as an input to a system which checks if the speed of vehicle is above the accepted level or not. For checking, the number in the speed limit image has to be extracted.

Can anyone suggest a simple method to perform number segmentation and identification in the above case?
Also if K-means algorithm can be used?

As with all ocr, you need to segment each individual character. Assuming this step is already done (and maybe some post processing to One easy way is to take the number and reduce the resolution. Similar to a lcd character display.

You can then do a simple correlation between the reduced resolution number from the image, and a template. Since the numbers on street signs are printed and very consistent a simple template match could work well.

your template should contain the numbers 0-9 in a font that closely matches those used on your road signs.


Learner level information on Kalman filtering for different input kinds I am learning Kalman filters and have seen example on data as state varaibles that have real values / numeric. However, in digital communication the information is in digital - bits. So, can Kalman filter not work on bits?

I have not found this information in any text book, tutorials and other reading materials. In text books and in wikipedia, it says that the state $x$ of the system is represented as a vector of real numbers I want to know if Kalman filter and Ricatti Equations needs to modified if the input signal $x$ is a discrete, taking values from a symbol set. Examples of such input is the DNA genes, BPSK symbols, quantized pixel values. Is there a different Kalman filter for symbolic input? Or can we use the same Kalman filter and its equations. Thank you

I am most convinced you should try Hidden Markov Chains (which i am not explaining here), just justifying its use:

Kalman Filter: continuous state space, discrete observations

https://stats.stackexchange.com/questions/183118/difference-between-hidden-markov-models-and-particle-filter-and-kalman-filter


Power of a speech signal I have been asked to calculate the power of a speech signal given by the following:
$$\ p(x)= \left(\frac a2\right)e^{-a\mid x \mid}$$

I know the formula for the power of a signal to be the following:

$$\frac{1}{T}\int_\frac{-T}{2}^\frac{T}{2}p(x)^2dx$$

Only this is not a periodic signal so I don't know where to set the upper and lower limits. I have tried taking the Fourier transform and evaluating but that gives me $\ P(f) = \frac{a^2}{a^2+2\pi f^2} $. If anyone would be able to guide me through a solution to this I'd really appreciate it. Simon.

for negative 'a' you can limit the integration up-to the 5-6 times of time constant of the dying exponential(as per your equation). this will give you the power approximately equal to actual.
The time constant of dying exponential is given by 1/abs(a). more details could be checked at https://en.wikipedia.org/wiki/Exponential_decay


Sampling - Higher order harmonics I have been studying about the sampling theorem and it seems that even though we sample at a frequency with the nyquist criterion, the harmonics (due to sampling process) remain within the nyquist frequency range, as illustrated below:

If we sample a 5MHz sine wave with 20 MHz sampling frequency (nyquist frequency being 10 MHz), its fundamental and harmonic frequencies (up to 4th) lie at:

fin=5 MHz (fundamental)

1st Harmonics: (around fs)
fs+fin=25 MHz
fs-fin=15 MHz

2nd Harmonics: (outside nyquist)
fs+(2*fin)=30 MHz
fs-(2*fin)=10 MHz

3rd Harmonics:
fs+(3*fin)=35 MHz
fs-(3*fin)=5 MHz (Inside nyquist)

4th Harmonics:
fs+(4*fin)=40 MHz
fs-(4*fin)=0 (DC)

I understand that when doing a high level simulation in matlab, we do windowing or coherent sampling to overcome this harmonics falling into the nyquist range or signal band by choosing a prime number of cycles.

But practically, if we give an arbitrary frequency to an ADC, which usually has a sample/hold circuitry in the front, aren't these harmonics generated by it?

Or do we have the control of the harmonics' power so that we can adjust it to a minimum (depending on the application).

(This question also relates to my previous question - Higher order harmonics during sampling)

I don't quite understand why you feel "harmonics" are relevant to this discussion.

A 5 Mhz has sine wave has no harmonics. If the signal has harmonics, it is not a sine wave any more but a different signal (rectangular, triangle, etc.)
For ANY signal: you need to determine the highest frequency that's in the signal and then chose the Nyquist frequency to be higher. Otherwise you get either loss of information or aliasing
Most ADC do have an anti-aliasing low-pass filter build in: if you have harmonics higher than Nyquist, the filter will just chop them off.

Points 2 and 3 apply to ANY signal: noise, music, video, OFDM, and also harmonic signals like rectangle our saw-tooth waves. There really is nothing special about "harmonics".


Bag of Tricks for Denoising Signals While Maintaining Sharp Transitions I know this is signal dependent, but when facing a new noisy signal what is your bag of tricks for trying to denoise a signal while maintaining sharp transitions (e.g. so any sort of simple averaging, i.e. convolving with a gaussian, is out). I often find myself facing this question and don't feel like I know what I ought to be trying (besides splines, but they can seriously knock down the right kind of sharp transition as well).

P.S. As a side note, if you know some good methods using wavelets, let me know what it is. Seems like they have a lot of potential in this area, but while there are some papers in the 90s with enough citations to suggest the paper's method turned out well, I can't find anything about what methods ended up winning out as top candidates in the intervening years. Surely some methods turned out to be generally "first things to try" since then.

L1 norm minimization (compressed sensing) can do a relative better job than conventional Fourier denoising in terms of preserving edges.

The procedure is to minimize an objective function

$$
    |x-y|^2 + b|f(y)|
$$

where $x$ is the noisy signal, $y$ is the denoised signal, $b$ is the regularziation parameter, and $|f(y)|$ is some L1 norm penalty. Denoising is accomplished by finding the solution $y$ to this optimization problem, and $b$ depends on the noise level.

To preserve edges, depending on the signal $y$, you can choose different penalties such that $f(y)$ is sparse (the spirit of compressed sensing):

if $y$ is piece-wise, $f(y)$ can be total variation (TV) penalty;
if $y$ is curve-like (e.g. Sinogram), $f(y)$ can be the expansion coefficients of $y$ with respect to curvelets. (This is for 2D/3D signals, not 1D);
if $y$ has isotropic singularities (edges), $f(y)$ can be the expansion coefficients of $y$ with respect to wavelets.

When $f(y)$ is the the expansion coefficients with respect to some basis functions (like curvelet/wavelet above), solving the optimization problem is equivalent to thresholding the expansion coefficients.

Note this approach can also be applied to deconvolution in which the objective function becomes to $|x-Hy|+ b|f(y)|$, where $H$ is the convolution operator.


Impulse response $h[n]$ from input and output signal using cross correlation on Matlab I have two signals on Matlab, one being input and the other output. I want to know if taking the cross correlation of these signals would give me the impulse response?And if so, will taking the Fourier transform of that resulting signal give me the frequency response?

i want to know if taking the cross correlation of these signals would give me the impulse response?

No.
Let's look at it in the frequency domain. Let's call the input $x$, the output $y$ and the impulse response $h$. Uppercase letters are spectra, lower case time domain signals.
For the transfer function we simply have $Y = H\cdot X, H = Y/X$. The spectra of the autocorrelation $r_{xy} = y \star x$ is
$$R_{xy} = X' \cdot Y = X' \cdot H \cdot X = H \cdot |X|^2$$
So you will only the get transfer function if the magnitude of the spectrum of the input signal is $1$ at all frequencies.
If you do time discrete auto correlation you also have to properly manage circular versus linear auto-correlation  .

and if so, will taking the Fourier transform of that resulting signal give me the frequency response?

The transfer function is indeed the Fourier transform  of the impulse response but since the cross correlation is NOT the impulse response, you won't be getting the transfer function either.


Advantages/disadvantage of zero padding I would like to know what are the advantages/disadvantage of zero padding with respect to frequency measurement and amplitude measurement

Advantages of zero padding:

If length of your sequence doesn't correspond to the size that can be handled efficiently with FFT routine (usually powers of prime numbers) then you might want to add some extra zeros to the nearest power in order to get the maximum speed-up. In worst case you double the memory you need for your signal.
Adding zeros is equal to interpolating samples of your spectrum with $\mathrm{sinc}$ function. Therefore you will find it looking more smoothly. It is not affecting frequency resolution in any way.
If some peaks were split between two bins, these can be interpolated to some extend and you might be able to retrieve some amplitude information thanks to that.
If you are trying to plot frequency response of some FIR filter from it's impulse response then you need to add zeros to it. Otherwise you will get as many samples as the length of filter is. Please see here
If you are trying to convolve your signal with some pattern using FFT, then you need to pad your signals with zeros to the appropriate length. Otherwise result of convolution will be incorrect (replicas in frequency domain will overlap). Please see here

Disadvantages of zero padding:

Can't think of any. Oh maybe except of people believing that it will improve their frequency resolution in some magical way - I knew PhD's saying that...


How to reduce synthesis artifacts produced by phase vocoder? I implemented the phase vocoder algorithm in Python to time-stretch speech signals by following this paper and referring to this MATLAB tutorial. I can distinguish words in the original signal from the resynthesis result, and the pitch, as expected, doesn't get altered even though the time axis has been stretched(which is the goal of using phase vocoder). However, there exist noticeable artifacts. Specifically, the output contains more noise resembling Gaussian noise compared to the input, and the words somewhat sound like being spoken a robot. I have tried some improvements such as identity phase locking suggested by this paper, but the result remained unchanged. Has anyone encountered similar issues? How can I improve the result?

The second artifact can be caused by the fact that phase vocoder slow-down resynthesis can reduce the frequency of any natural tremolo and vibrato in the human speaking voice, as well as any tiny period-to-period jitter in glottal closures.

One possible solution might be to measure the frequencies of tremolo and vibrato, and the jitter statistics, of voiced windows, and re-add those artifacts to (or modulate) the result of the PV resynthesis, if they are absent or too small.


General Questions MIMO Communication theory The following are questions regarding MIMO communication theory. I have come across a MATLAB implementation of MIMO channel. My questions are not regarding the implementation, but rather regarding the theory itself.
In MATLAB there is comm.MIMOChannel that creates a MIMO channel. The following are the properties which one can play with

 MIMOChannel properties:
%
%   SampleRate                  - Input signal sample rate (Hz)
%   PathDelays                  - Discrete path delay vector (s)
%   AveragePathGains            - Average path gain vector (dB)
%   NormalizePathGains          - Normalize path gains (logical)
%   FadingDistribution          - Rayleigh or Rician fading
%   KFactor                     - Rician K-factor scalar or vector (linear scale)
%   DirectPathDopplerShift      - Doppler shift(s) of line-of-sight component(s) (Hz)
%   DirectPathInitialPhase      - Initial phase(s) of line-of-sight component(s) (rad)
%   MaximumDopplerShift         - Maximum Doppler shift (Hz)
%   DopplerSpectrum             - Doppler spectrum
%   SpatialCorrelation          - Spatially correlated
%   NumTransmitAntennas         - Number of transmit antennas
%   NumReceiveAntennas          - Number of receive antennas
%   TransmitCorrelationMatrix   - Transmit correlation matrix (or 3-D array)
%   ReceiveCorrelationMatrix    - Receive correlation matrix (or 3-D array)
%   AntennaSelection            - Optional transmit and/or receive antenna selection
%   NormalizeChannelOutputs     - Normalize channel outputs (logical)
%   RandomStream                - Source of random number stream
%   Seed                        - Initial seed of mt19937ar random number stream
%   PathGainsOutputPort         - Enable path gain output (logical)

I have the following questions: ( I really hope I can get some help)

Does the sample rate mean the actual symbol rate of the input? Is it related to the modulation we use?
Do path delays dictate whether we have a frequency selective or frequency flat fading channel?
The channel above gives the option of having either Rayleigh or Ricean fading channel gains, I know that Rayleigh is good approximation because of the central limit theorem (i.e large number of random paths), if one specifies that the distribution is Rayleigh, does that mean that the Path Delays should be infinite? Or does this distribution mean a distribution of the path itself.
Why would a Direct PATH (LOS) have a doppler shift different that the maximum doppler shift? Then which one dictates the coherence time?

Thank you for any advice or answers on any of the above questions

As it says in the documentation, SampleRate is the signal's sampling rate, not the symbol rate.
The path delays together with your symbol rate determine whether the channel is flat or frequency-selective.
You should choose Rayleigh if there is no line-of-sight component in the received signal, and Rice otherwise. Regarding the path delays: the model is that the gain you observe at each sampling instant is the result of many paths, whose arrival times are so close that you can't resolve them. So, in your model, you always have a finite number of paths, each of which has random Gaussian gain.
Each path has a different length and arrival angle, so it can have a different Doppler.

As a side note, and expressing only my opinion, it is more educational to try to create your own models of the fading channel than using a Matlab toolbox. Toolboxes are created to increase an expert's productivity, not for learning.


FMCW Radar - Range-Doppler Plot The FMCW radar with the following parameters was implemented:
% System parameters Value

% ----------------------------------

% Operating frequency (GHz) 77

% Maximum target range (m) 15

% Range resolution (m) 0.037

% Maximum target speed (km/h) 30 for fall detection

% Sweep time (microseconds) 50

% Sweep bandwidth (GHz) 4

% Maximum beat frequency (MHz) 73

% Sample rate (GHz) 4.050

A target was simulated with the following parameters(static radar)

dist = 14;

speed = 10*1000/3600; ;%2.77km/h

rcs = 100;

and the Range Doppler image plotted using the phased array functions.

RangeFFTLength = 2048;

DopplerFFTLength = 512;

rngdopresp = phased.RangeDopplerResponse('PropagationSpeed',c,... 'DopplerOutput','Speed','OperatingFrequency',fc,'SampleRate',fs,... 'RangeMethod','FFT','SweepSlope',sweep_slope,... 'RangeFFTLengthSource','Property','RangeFFTLength',RangeFFTLength,... 'DopplerFFTLengthSource','Property','DopplerFFTLength',DopplerFFTLength);

figure

plotResponse(rngdopresp,xr);

[![enter image description here][1]][1]

The results are as expected.

Interested to get the same results implementing the 2DFFT manually.

For this purpose the following code was used:

Nr = RangeFFTLength;

Nd = DopplerFFTLength;

% 2D FFT using the FFT size for both dimensions.

sig_fft2 = fft2(xr,Nr,Nd);

% Taking just one side of signal from Range dimension.

sig_fft2 = sig_fft2(1:Nr/2,1:Nd);

sig_fft2 = fftshift (sig_fft2);

RDM = abs(sig_fft2);

RDM = 10*log10(RDM);

Yielding the following:

The index of the range and velocity are: [518 229].
This yields the following values:
Range:
fb=fs/2;

range_max = cfb/(2sweep_slope);

distance = ((518-512)/512)*range_max

distance =

43.9453

This is 3 times larger then the correct value. Why?
Velocity:
Converting the index to the abs value using the following, yields correct value:
vel_max = lambda/(4*tm);

sp = ((293-256)/256)*vel_max

sp =

2.8155

Seems, however, that the velocity in this case is positive – why?
The code can be seen below:
https://drive.google.com/file/d/1e7E6ItgVeHTaU1qx3VL_aDfB7nes7ipq/view?usp=sharing

So the problem is that when simulating FMCW signals, there is some ambiguity in the convention on the mixing order. MATLAB's documentation makes reference to it here:

"The mix operation reverses the Doppler shift embedded in x, because
of the mixing order of xref and x. The mixing order affects the sign
of the imaginary part of the output argument, y. There is no
consistent convention in the literature about the mixing order. This
function and the beat2range function use the same convention. If your
program processes the output of dechirp in other ways, take the mixing
order into account."

So you're amost there. MATLAB's functions account for this so you need to account for the change yourself when doing the manual range-Doppler map. You can do this by applying a complex exponential to each pulse that shifts the phase in the other direction so that it will produce the right Doppler shift.
For example, in the simulation loop right before you assign the dechirped signal to xr do something like the code below. I currently don't have the Phased Array Toolbox installed so I did some of the calculations manually:
for m = 1:Nsweep

    % After all the target sim stuff do your own Doppler calculations
    rangeToTarget = norm(tgt_pos - radar_pos);
    radialVel = dot(tgt_vel - radar_vel, tgt_pos - radar_pos)/rangeToTarget;

    % Negate Doppler so that closing targets produce positive frequencies
    fd = -2*radialVel/lambda;

    % De-chirp and apply phase shift for this pulse to reverse Doppler
    % to the correct value
    xr(:,m) = dechirpsig*exp(1i*2*pi*(2*fd)*(m - 1)*tm);

end

This uses the simple "stop and hop" model where the target moves in discrete time steps which are equal to the sweep time tm.
Addressing Range Issue
I totally to forgot to address the range part of your question. The reason this is happening is because when you use MATLAB's functions to generate the range-Doppler map, you are giving it a sampling rate $f_s = 4 \space GHz$, which in turn yields a range-bin (or range-cell) size $\delta R $ of
$$\delta R = \frac{c}{2B} = \frac{c}{2(4e9)} = 0.0375 \space m$$
If you had set RangeFFTLength = numel(waveform) when doing your own 2D FFT, you would have a HUGE number of samples in the range dimension, but it would be correct given the achieved range-bin resolution $\delta R$ and will give you the correct target ranges.
The problem now comes that when you use less range samples than what you started with by setting RangeFFTLength = 2048. This decimates the effective sampling rate, so your new range-bin size is now going to be artificially larger. If you take the ratio of the number of samples produced by using $f_s = 4 \space GHz$ and the size of the range FFT of 2048, your effective sampling rate is now 98 times worse! So your new range bin size becomes
$$\delta R_{new} = 98\frac{c}{2B} = 98\frac{c}{2(4e9)} = 3.675 \space m$$
Applying this range-bin size to the peak range index should yield the correct range in your range-Doppler map. However, note that this is still wrong. I say that because the whole point of putting that much bandwidth into a chirp is so that each range-bin in the range-Doppler map is equal to to 0.0375 m and NOT 3.675 m.
If you did it correctly, you would be able to sample at a much lower rate since your beat frequencies after mixing have lower bandwidth requirements, giving you a smaller and reasonable number of samples. Each range-bin size is still 0.0375 m, preserving the range resolution desired.


How to correctly classify voice from other sound events? I have a huge labeled dataset of several thousand sound events, including human voice, dish washing, things falling to the ground, among others.

I need to report when a human voice event takes place. Notice that it's not a simple voice activity detection (VAD), since there are other sound events competing with voice (in VAD, they normally only concern about well behaved noise).

My approach so far has been to train a binary SVM classifier (voice and non-voice classes) with MFCC features. Even after parameter optimization and tinkering about with different number of MFCC coefficients, the performance is awful for such a simple task...

Is there any heuristic or any thing that may help distinguish voice from non-voice events that I'm missing out on?

(This related question is similar, but I don't need to "eliminate" other sounds completely, and I'm looking, first, for possible heuristics to improve the classification. This article would be my last resort.)

It is very possible that you don't have enough data to train your classifier.  Although several thousand seems like a lot, it means that you only have a fraction of the actual thing you are trying to identify.  It may be helpful to expand your training data by slightly tweaking your samples to improve the robustness of your machine learning algorithm.  You could take the existing human voice samples you have and try slowing each one down slightly, increasing the speed slightly, introducing pink/white noise, etc...  This technique is commonly used with images in CNNs (by introducing small rotations and stretching different ways) to increase the size of the training set when training image classifiers.

Also, don't forget to reduce your learning rate as learning proceeds.  With too large of a step size your gradient decent could actually overshoot the loss function minima and result in poor results or saturate learning prematurely.

Hope that is helpful.


Algorithm to identify pitch-shifted recordings? new to this whole DSP thing. Two related questions:

Does there exist an algorithm or classifier that, when given two sound files, would be able to tell me if one is a pitch-shifted version of the other? There are algorithms that can tell you whether two images are resized versions of one another, I feel like this would be the audio equivalent.
How would you write an algorithm that takes two recordings that are pitch-shifted versions of one another and spits out the amount of pitch shift. My naive guess would be to count zero-crossings in a fixed window for both recordings and compare the two counts, but I feel like that might lead to aliasing issues.

Thanks!

There is an algorithm called Dynamic Time Warping that allows you to measure a distance between 2 time-series indepentently of their relative time stretching. Thus with a correct threshold value, you could detect if the two recording (the original and the pitch-shifted) are actually the same.

Then, by using the warping path, you could determine the relative amount of pitch shift between the 2 recordings.


What's the relation between frequency band of $X(j\omega)$ and $\Phi_{xx}(j\omega)$? in which:

$x_{c}(t)$ is a continuous-time signal
$X(j\Omega)$ is the Fourier Transform of $x_{c}(t)$
$\Phi_{xx}(j\Omega)$ is the Power Spectrum Density of $x_{c}(t)$ which defined as Fourier Transform of auto-correlation of $x_{c}(t)$.

in other words I want to know when it is said that a signal is band-limited which is the case? band-limited according to $X(j\Omega)$ or $\Phi_{xx}(j\Omega)$?
or if it's a relationship between two cases what it is.
thanks.

The autocorrelation of $x(t)$ is
$$r_x(t)=x(t)\star x(-t)\tag{1}$$
where $\star$ denotes convolution. Taking the Fourier transform of $(1)$ gives
$$S_x(\omega)=X(\omega)X^*(\omega)=|X(\omega)|^2\tag{2}$$
$S_x(\omega)$ is the energy density of $x(t)$, and according to $(2)$ it equals the squared magnitude of the Fourier transform of $x(t)$. So if $x(t)$ is band-limited, both $X(\omega)$ and $S_x(\omega)$ are zero outside the signal's bandwidth.
Note that a deterministic continuous signal which has a Fourier transform (represented by an ordinary function) is usually an energy signal, which doesn't have a power spectrum (only an energy density).


Can software distinguish between different musical instruments If possible, I would like to write a little script that can identify different musical instruments from just a recording.  Apart from variations of amplitude and frequency over time (which can be shown by time-dependant and spectrogram, respectively), what other variables are there?

Furthermore, how would one do the classification?  Fingerprinting is an option, but that would only work if a single piece at a time is played.

What you want to do is timbre recognition.


How can I increase image size by zero padding? For increasing the image size from 8X8 to 32X32, I tried to do zero padding. The simplest concept of zero padding in Java for me seems to be the code as shown below in Java such that I can shift each point by 4 points while rescaling.

for (int i = 0; i < imageWidth; i += 4)
        for (int j = 0; j < imageHeight; j += 4)
            pspectrum[i][j] = powerspectrum[(i + 1)/4][(j + 1)/4];

Now, what I am mainly concerned is that is there any more accurate method and how this helps with my calculation.

I suppose the array names are not descriptive of the contents and you have simply repurposed them. To zero-pad from 8x8 to 32x32:

for (int i = 0; i < 32; i++)
    for (int j = 0; j < 32; j++)
        pspectrum[i][j] = 0;
for (int i = 0; i < 8; i++)
    for (int j = 0; j < 8; j++)
        pspectrum[i][j] = powerspectrum[i][j];

The purpose of this is that the discrete Fourier transform (DFT) of a larger image has a higher frequency resolution, so you can resolve for example the dominant frequency more accurately. It is the same as really well done interpolation in the frequency domain.


iZotope RX's FFT processing? How does iZotope RX handle the FFT (i.e. how does it manage to perform those processes in the frequency domain without creating glitches) since its able to e.g.

Apply another DSP process (even time-domain) into a spectral part (i.e. apply e.g. a phaser onto individual parts of the FFT).
Attenuate or gain spectral    parts (without touching the surrounding parts).

Is it really frequency domain filtering or time domain using time domain filters?

I thought this kind of processing was impossible: Why is it a bad idea to filter by zeroing out FFT bins?

Spectral Repair is more complex than just muting a handful of harmonics.

But first, let's just clarify that as far as Linear Time Invariant (LTI) systems are concerned, their processing can be applied, equivalently, either in the time domain (via the operation of convolution) or the frequency domain via either the overlap-add or overlap-save methods.

iZotope's Spectral Repair performs interpolation. That is, to remove an unwanted portion of a recording, it will try to "mimic" the characteristics of the spectrum in either ends of the unwanted region.

Consider the following image:

Can't we just remove the railings by copying and pasting the grass from another (and very similar) part of the image and adjust its brightness? Couldn't we even replicate how the green colour varies in the grass portions and just generate some "new" grass for the regions of the railings?

Is it going to look odd? Possibly, if it was overdone, but for small regions it may be good enough to fool the eye.

This is what the "Replace" mode (and "Attenuate") does in iZotope, but instead of grass and railing there is background noise and local disturbance.

It's not so much attenuation (or "setting harmonics to zero") as "masking" or "hiding" the unwanted sound, since it attempts to make up a good enough patch of harmonics to bury it in by "looking at" the surroundings of the disturbance.

For more information on interpolation please see this link. For an example of "learning the profile and applying a filter" please see this link and this link.

Hope this helps.


Is it possible to use a control scheme with a delay as prefilter? I am studying control systems, and in my studies I have encountered the topic of control scheme of the form

and I have seen thaht this control scheme is really useful since it is possible to add a lag or lead compensator outside the closed loop, so it won't give stability problems.

At this point I was curious to know a little more and in particular I would like to know if it is possible to use a prefilter which  is a delay. In such a way it would be possible to avoid the problem of stability given by a system with delay.

So, to explain better, I would have for example:

$T(s)=e^{-s}\frac{C(s)P(s)}{1+C(s)P(s)}$

and so I could design a stable closed loop without delay, and then filter the reference with a delay, so use:

$C_r(s)=e^{-s}$

instead of using a closed loop with delay, which gives stability problems.

Do you think it would be a good idea?  Or it would be impossible or conceptually wrong?

1 - It can be a good idea to put a filter $C_r(s)$ in order to smooth the reference. For example, if you have a speed regulator, and the reference jumps from 100 km/h to 200 km/h, having a low-pass filter will smooth the actual reference fed to the closed-loop. It could prevent excitation of unmodeled dynamics, actuator wind-up, etc. That being said if $C_r(s)$ is a simple delay, it has no purposes as far as I can see.

2 - $C_r(s)$ only affects the reference fed to the closed-loop. It will not affect stability of the closed-loop like you mentioned. It also has no effect on disturbance rejection or errors due to modelization (such as unmodeled dynamics).


Given the input and the output, how to determine the impulse response? I would like to find the impulse response, $h[n]$, of an LTI system given the input

$$x[n] = [1,-3,2]$$

and the output

$$y[n] = [1,-1,-4,4]$$

I know that $y[t]=x[t]*h[t]$, but I am having hard time to figure out the right way to calculate the impulse response. I know very little about signal processing, so if you don't mind giving an easy explanation, then I appreciate it. Or, if it's possible to do an example, that's better.

$\mathcal Z$-transforming the input and the output and then dividing them, we obtain the transfer function

$$\begin{array}{rl} H (z) &= \dfrac{1 - z^{-1} - 4 z^{-2} + 4 z^{-3}}{1 - 3 z^{-1} + 2 z^{-2}}\\\\ &= \dfrac{z^3 - z^2 - 4 z + 4}{z (z^2 - 3 z + 2)}\\\\ &= \dfrac{(z-1)  (z-2) (z+2)}{z (z-1) (z-2)}\\\\ &= \dfrac{z+2}{z} = \color{blue}{1 + 2 z^{-1}}\end{array}$$

Computing the inverse $\mathcal Z$-transform of $H (z)$, we finally obtain the impulse response.


Combining different Likelihoods - Particle Filter I am new to particle filters and have a problem when it comes to the weightening (using SIR to be precise).
The problem is that my different state-variables are subject to different noise-realizations (e.g. one variable experiences additive gaussian noise while another one is log-normally distributed), how would one normally deal with these situations?

What I tried out first was to calculate the probability of each state-variable of a particle separatly according to the respective likelihood, then multiply all resulting occurence proabilities for the respective particle and afterwards noramlize the resulting final proabilities as weights.

Sadly, this does not work very well (I suspect this might be the case because it ignores possible correlations between the state-variables) so I am very thankful for any kind of advise in this matter.

EDIT:
From a more mathematical perspective, what I do is
$$f(\mathbf{z}\mid \mathbf{x})\approx f(z_{1},z_{2}\mid \mathbf{x})=f(z_{1}\mid \mathbf{x})f(z_{2}\mid z_{1},\mathbf{x})\approx f(z_{1}\mid \mathbf{x})f(z_{2}\mid \mathbf{x})$$

So first I assume that the multivariate distribution can be expressed by a distribution of its components, then in the second approximation step I ignore the correlation between the respective components. So I would like to know under which restrictions the first step is allowed (I for example know  that it its correct for a multivariate normal distribution) and also how I could enhance the second approximation $f(z_{1}\mid \mathbf{x})f(z_{2}\mid \mathbf{x})$ to balance the lack of the correlation.

Let $\mathbf z = [z_1, z_2]^T$ be the observations and $\mathbf x = [x_1, x_2]^T$ be the hidden states, and the observation model is $z_i = x_i + n_i$ for $i=1,2$. Here $n_i$'s are independent (but not identical, since the noise has different distributions).

  So I would like to know under which restrictions the first step is allowed (I for example know that it its correct for a multivariate normal distribution)

By definition, the density function of a random vector is given by the joint density of its components. So $f(\mathbf z | \mathbf x) = f(z_1, z_2 | \mathbf x)$ is always true.

  and also how I could enhance the second approximation

The second approximation is actually an equality under the observation model with independent components. Since $z_2 = x_2 + n_2$, conditioned on $x_2$, we have $z_2 \perp (z_1, x_1)$. This implies that $f(z_2 | \mathbf x, z_1) = f(z_2 | x_1, x_2, z_1) = f(z_2|x_2)$. And similarly, $f(z_1|\mathbf x) = f(z_1 | x_1)$.


Why do we use the HSV colour space so often in vision and image processing? I see the HSV colour space used all over the place: for tracking, human detection, etc...  I'm wondering, why?  What is it about this colour space that makes it better than using RGB?

The simple answer is that unlike RGB, HSV separates luma, or the image intensity, from chroma or the color information.  This is very useful in many applications.  For example, if you want to do histogram equalization of a color image, you probably want to do that only on the intensity component, and leave the color components alone.  Otherwise you will get very strange colors.

In computer vision you often want to separate color components from intensity for various reasons, such as robustness to lighting changes, or removing shadows.

Note, however, that HSV is one of many color spaces that separate color from intensity (See YCbCr, Lab, etc.).  HSV is often used simply because the code for converting between RGB and HSV is widely available and can also be easily implemented.  For example, the Image Processing Toolbox for MATLAB includes functions rgb2hsv and hsv2rgb.


How to predict the cramped frequency of a digital filter based on an analogue frequency? I have a second order analogue high pass transfer function (unity gain at infinity). It's magnitude response hits the -20 decibel line at a frequency of 5706 Hz (the corner frequency is 18000 and sample rate is 44100). When I convert this analogue filter to a digital IIR via the BLT method, the digital filter's magnitude response hits the -20 line at 11149 Hz. This will happen with the standard RBJ audio EQ cookbook formula i.e. it's not a problem with my maths. I know that digital filters suffer from cramping and that's what I'm interested in. Is there a formula to predict what an analogue frequency will be cramped to when the digital filter is realized given the analogue filter's frequency? i.e. something that will tell me 11149 given 5706. Just to be super clear I'm not asking how to do pre-warping using: $ tan(\frac{w_c}{2}) $.

I've come up with a solution.

Let $ f_a = $ analogue frequency, $ f_d $ = digital frequency, $ f_s = $ sampling rate and $ f_c = $ corner frequency with $ \omega_c = \frac{\pi f_c}{f_s} $ and $ \omega_a = \frac{\pi f_a}{f_s} $ then:

$$ c = \frac{\omega_c}{\tan(\omega_c)} $$

$$ \omega_d = \arctan\left(\frac{\omega_a}{c}\right) $$

Giving:

$$ f_d = \frac{\omega_d f_s}{\pi} $$


What can be a limit to the maximum possible power of a pilot signal In practice, what can be a limit to the maximum power of a pilot signal in wireless communication?

There is no limit other than what is imposed by the FCC and underlying hardware. Any commercial/consumer device has an operating power defined either on the device or in the manual. The max power of the signal (and pilot) will be within those bounds.


Image processing books Gonzalez et al. vs Sonka et al? I don't know where else to post this question so here it goes. I'm trying to buy a DIP book and the best two right now are:

Image Processing, Analysis, and Machine Vision by Roger Boyle, Milan Sonka, Vaclav Hlavac 4th ed.
Digital Image Processing by Rafael C. Gonzalez, Richard E. Woods 4th ed.

They're both pricey books. So if anyone has read these books and can make a recommendation or a pros/cons comparisons of the books, I'd appreciate it.

I've already bought the "Digital Signal Processing: A Practical Guide for Engineers and Scientists by Steven W. Smith".

I am not sure there are "best" image processing or computer vision books to buy. The topics covered can be very wide, and some can be better on some aspect (morphology, segmentation, denoising, etc.).

Some sites recommend a hanful of such books, and link to their draft versions, as in 8 Books for Getting Started With Computer Vision. Indeed, a lot of authors have a draft on a personal page. This is the case for instance for Computer Vision: Algorithms and Applications by Richard Szeliski.

CVonline: Vision Related Books including Online Books and Book Support Sites also points to classics, like:

D. H. Ballard, C. M. Brown; Computer Vision, Prentice-Hall Inc New Jersey, 1982, ISBN 0-13-165316-4.
A. Blake, M. Isard; Active contours, Springer, London, 1998, ISBN 3540762175.


Linearized system and State Space I want to ask about this scheme:

$u_0$ should be something as input in specific point given by the initial conditions and $y_0$ output. Block State Spase represented linearized system by A, B, C, D state metrices.

And I'm wondering why there is "-" by $u_0$ and "+" by $y_0$. Is this scheme universal? Why does it look like it looks?

The scheme is not universal.

One reason to do it is to allow the system to be modeled around a set-point or bias level ($u_0$) without having that set-point affect the state space model.

Subtracting the set-point, doing the state space equations, and then adding back on the corrected set-point ($y_0$) lets you do this.


How to change the FFT phase domain of each frame by 90 degrees I have read a signal in Matlab, resampled and implemented an FFT of each frame on it, now I want to rotate the domain of phase in each frame 90 degrees.

I have a vector of FFT and used the function named rot90(matrix), that rotates a matrix 90 degrees, but I'm not sure it is what I wanted to work correctly on my function.

How can I do that in some other way?  Or does the rot90 correctly work on changing the phase of domain in FFT which has been imported into a vector?

It sounds like you're trying to plot the FFT such that the result is centered at 0, instead of at Fs. Matlab has a function for this fftshift.

Here's a simple example when trying to plot either the onesided or twosided FFT.

f = 100; % Signal frequency (Hz)
fs = 1024; % Sampling frequency (Hz)
N = 1024; % Number of samples
Nfft = 1024; % FFT size

xf = 20*log10(abs(fft(x, Nfft))); % Compute the PSD

fw1 = fs * (0:Nfft/2)./Nfft; % Onesided domain
fw2 = fs * (-Nfft/2:Nfft/2-1)./Nfft; % Twosided domain

figure
subplot(211), plot(fw1, xf(1:Nfft/2+1)), title('Onesided')
subplot(212), plot(fw2, fftshift(xf)), title('Twosided')


Direct Inverse Filtering an image, spatially convolved with a blurring function without any noise, doees not give my original image back? I am trying to learn about inverse filtering of images and am trying it out on a test image in MATLAB.

The test image is an 8x8 checkerboard, with 64x64 pixels, generated by the MATLAB code:

g = checkerboard(8);

I am using a 3x3 averaging filter, h, and convolving it with the test image.

I am not adding noise and I want to see the results of direct inverse filtering the distorted image.

I am applying the deconvolution function in the frequency domain by:

1) Padding the 3x3 averaging mask h to generate a 64x64 mask hp.

2) Centering the coefficients (0.111...) of hp to maintain symmetry.

3) Multiplying hp by (-1)^x+y (since I want the fft of the mask to give me a centered Fourier transform, hpc.

4) Taking the FFt of the padded, centered mask and multiplying (-1)^(u+v) (since I centered the mask in the spatial domain in step 2) to give me H.

5) Taking the distorted, and frequency shifted test image and dividing it by H, (ie. Y = G./H). (Note: I am not padding the image since the mask was padded, and wraparound is not a huge concern here)

6) Obtaining the inverse filtered image inv = (-1)^x+y * real(ifft2(Y))

My question: Why can't I obtain the original (or something very close to it)?

I think I found out the answer to my own question.

Since I'm using an averaging filter mask (some constant multiplied by a matrix of ones), the Fourier Transform of such a filter is an impulse at zero frequency.

When I pad the mask filter in step 1 to make the dimensions of the mask equal to the image, I am 'truncating' the filter, ie. windowing the mask, which is akin to multiplying the Fourier Transform of the mask with a sinc function.

As a result, the Frequency Domain filter I use in step 5, ie. H, is not an impulse but a sinc function. That's the reason why my results are not ideal.

On the other hand, if I were to use an impulse in step 5, only the zero-frequency value would remain, and my result would be a constant image.

In short, it's a lose-lose situation.