Information to be transmitted being a stochastic process? I am reading Digital Communication, 3rd edition by Barry & Lee,

Book's Webpage:
enter link description here

This is an excerpt taken from this book (mentioned above) from Chapter 3, Page # 57, 1st Para, 1st Line.

It says,

 Although modulation and demodulation are deterministic,the information to be
 transmitted over a communication system, as well as the noise encountered in the
 physical transmission medium,is random or stochastic.

I understand Noise being a random process but how the information to be transmitted could be a random process.
I mean do we not know the nature of information to be transmitted by us ahead of any transmission ?

Please elucidate this.

The excerpt is intimating that, for communication applications, the information source is stochastic in the sense that its output is not fully known a priori. If you did know its output, then you wouldn't need a communication link at all! With that said, there are levels of gray to this effect, and there is a large area of study on its own known as information theory. Here's a very basic high-level summary:

The amount of information conveyed by a transmission is related to the amount of entropy contained in the transmission. Its entropy measures how unpredictable the transmission is. In your example, if a receiver observes a signal that is fully deterministic (i.e. it is completely predictable), then that transmission contains zero information. So, in a sense, the less we know what to expect about the next observed transmission, the more information that it can carry.

So, in terms of maximizing information rate, you want to maximize the entropy in your transmitted signal. However, there are some cases where reducing your transmission's entropy could be advantageous. Specifically, error-correcting codes are encoding methods that allow a receiver to recover from symbol errors by exploiting intentionally-inserted redundancy into a transmitted signal (which reduces its entropy). Thus, the information rate is lower, but the probability of error is reduced as well.


Solving Convex Optimization Problem Used for High Quality Denoising The highest voted answer to this question suggests that to denoise a signal while preserving sharp transitions one should

  minimize the 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.

However, there is no indication of how one might accomplish that in practice as that is a problem in a very high dimensional space, especially if the signal is e.g. 10 million samples long. In practice how is this sort of problem solved computationally for large signals?

Boyd has A Matlab Solver for Large-Scale ℓ1-Regularized Least Squares Problems. The problem formulation in there is slightly different, but the method can be applied for the problem.

Classical majorization-minimization approach also works well. This corresponds to iteratively perform soft-thresholding (for TV, clipping).

The solutions can be seen from the links. However, there are many methods to minimize these functionals by the extensive use of optimisation literature.

PS: As mentioned in other comments, FISTA will work well. Another 'very fast' algorithm family is primal-dual algorithms. You can see the interesting paper of Chambolle for an example, however there are plethora of research papers on primal-dual methods for linear inverse problem formulations.


Basics of waveform, how to define the zero point Audio is a sequence of waveforms which have peaks, troughs, 0, just like below.

How to define the x axis (y = 0) in the picture, how to choose the 0 point in the waveform since there are so many points on the waveform?

Why not move the waveform above or below the x axis(y=0)?

This question seems a little stupid. But I still hope someone can help with it.

You can't produce a DC component with speech, so the line $y=0$ is defined by the speech signal having zero mean.


IIR oscillator produces stable output in MATLAB, unstable on chip, how do I tune it? I'm trying to make an oscillator using IIR filter and I'm having problems tuning it so that it's stable on the platform I'm using (ARM Cortex M4F and CMSIS DSPLib's arm_biquad_cascade_df1_f32). Note that this is my first time doing this and that I don't know much about control theory, so I could be missing something very obvious here.

Looking through some textbooks, I found the following formula for a transfer function of an oscillator:

$$H(z)=\frac{Asin(\omega_0)}{1-(2rcos(\omega_0))z^{-1}+r^2z^{-2}}  $$
where, as far as I understand it, A is the amplitude of the oscillator, $\omega_0$ is the normalized angular frequency, which is calculated as $\omega_0=\frac{2f_{real}}{f_{sampling}}$ and the r is parameter which sets the poles distance to the unit circle and decides if the system's output will be dampened, opposite of dampened (can't remember the word at the moment) or constant.

I designed an oscillator using following parameters I got from some other requirements:
sampling frequency $F_s=759493.7$, oscillator frequency $f1=25630$, amplitude $A=1$ (I plan to fine-tune this later), $r=1$.

This gives me then $\omega_0=\frac{2f_1}{F_s}=0.0674923$ and the transfer function of:
$$H(z)=\frac{0.0674411}{1-1.9954464z^{-1}+1z^{-2}} $$

Looking at the impulse response in MATLAB, I get the expected result:

So far so good. I then tried to implement this using arm_biquad_cascade_df1_f32. Here's the block diagram of the filter according to documentation:

Since I only have a single second order section in my oscillator, I'm using only one section as shown on the picture.
I set the coefficients to be:
$b_0=0.0674923$, $b_1=0$, $b_2=0$
$a_1=-1.9954464$, $a_2=1$

When I start executing this filter, the output is unstable and quickly goes to infinity. Manually experimenting with the parameter $r$, I managed to get a slowly decaying response with $r=0.5$. I suspect that what I'm seeing are the results of the MCU being only able to work with 32bit floating point numbers and lacking precision.

So my question is: How do I model the oscillator while taking into account the numerical precision errors?

A first thing: I think there is a $\pi$ missing in your conversion between frequencies and angular frequency. Your impulse response shows a period of about 90, while your period should be about 29 samples - there's clearly something wrong here.

What you are seeing on your target hardware is due to the difference in precision between matlab's implicit use of double precision (64-bit), and the use of 32-bit for the routines you chose to use.

Even if this is not very apparent through the transfer function or the biquad structure, what this "filter" with no input does is simply rotate its initial conditions. With 32-bit precision, there is no guarantee that the magnitude of the input vector will be preserved by rotation - it can either grow (instability) or decay.

If you intend to use this as an oscillator (that is to say, if the input will always be null except for the initial conditions):

Regularly - say every 100th sample - mess with your state variables to readjust the norm of the vector to 1.
Or use another simpler and more stable technique like a phase accumulator and a sine lookup table.


Understanding Eigenface Face Recognition Algorithm Notations I have been reading the paper from Matthew A. Tusk and Alex Pentland

https://www.cs.ucsb.edu/~mturk/Papers/mturk-CVPR91.pdf

I got stuck in this part. I have a couple of questions.

Assume my images are 256 by 256.

1- What are G1,G2...Gm ? Row vectors? Column vectors? (1*65536) or (65536 * 1) ? Or are they 2D vectors (256*256)?

2- Then it says A = [ F1 F2 ... Fm] . Now I do not get this notation. Again A is what ? A matrix of matrices? or just a huge column/row vector?

Thanks

G1, G2... Gm are column vectors of length 65536.  The matrix A has columns that are the principal components.


Root Locus for a system I know this subject tends more to control theory but I am certain its part of the global knowledge basis for engineers, thus I believe I can find the answer here.
I'm trying to draw Root Locus for $K<0$, but can't understand how the system behaves.

I have this open loop system:

$$
\frac{20(1+s)^2}{(1+\frac{s}{10})(1+\frac{s}{2})}
$$

I know it has no asymptotes, but yet by one of the drawing rules for $k<0$, it tends to $\pm \infty$ in the real axis. How is it possible?

I will denote the given system as,

$$
G(s) = \frac{20(1+s)^2}{(1+\frac{s}{2})(1+\frac{s}{10})}
$$

The input or reference signal to the closed loop will be denoted by $R(s)$ and the output by $Y(s)$. The controller will initially be denoted by $H(s)$, while in this case it is equal to $K$. The  block diagram of the closed loop system can then be drawn as,

With its transfer function from input to output equal to,

$$
\frac{Y(s)}{R(s)} = \frac{G(s)}{1 + H(s) G(s)}.
$$

For root locus analysis you want to find the poles of this transfer function as function of $H(s)$. The poles can be found by solving when the denominator is equal to zero. However the denominator (and numerator) are polynomials of $s$, so can't be fractions themselves. In order to ensure this I will slit $G(s)$ up in to two parts, its numerator, $G_n(s)$, and denominator, $G_d(s)$. The same can be done for $H(s)$. The resulting transfer function can then be written as,

$$
\frac{Y(s)}{R(s)} = \frac{\frac{G_n(s)}{G_d(s)}}{1 + \frac{H_n(s)}{H_d(s)} \frac{G_n(s)}{G_d(s)}} = \frac{H_d(s) G_n(s)}{H_d(s) G_d(s) + H_n(s) G_n(s)}.
$$

As stated before, in this case $H(s)$ is equal to $K$, thus $H_d(s)$ is equal to one. By substituting in $K$ for $H_n(s)$ and the two polynomials, from this first equation, for $G_n(s)$ and $G_d(s)$, then the following expression can be obtained,

$$
\frac{Y(s)}{R(s)} = \frac{20(1+s)^2}{(1+\frac{s}{2})(1+\frac{s}{10}) + 20K(1+s)^2}.
$$

The resulting denominator is a second order polynomial, so solving for the poles yields two solutions, which can be simplified to,

$$
p=-\frac{400 K \pm 4 \sqrt{1 - 225 K} + 6}{400 K + 1}.
$$

For $K=-\frac{1}{400}$ the limits of both $p$ are $\pm\infty$ depending if you approach that value from the left or the right. Thus near $K=-\frac{1}{400}$ there will be values for the poles of the closed loop transfer function which will have a positive real values.


Extrapolation of sinusoidal signal Let there be a discrete signal that is the sum of sinusoidals and can be described by
$$s(n)=\sum_i A_i \mathrm{sin}(2\pi f_i \frac{n}{f_s}+\phi_i)$$
where $A_i, f_i, \phi_i$ are unknown but fixed parameters of my signal, $f_s$ is my sampling rate, and $n$ is the sample measured at time $t=\frac{n}{f_s}$.
Of this signal I measure the first $N$ samples ($n=1,2,3...N$) and I also measure another $N$ samples at points $n=N+20000,....2N+20000$.
I am now interested in estimating/predicting/extrapolating the value of a single sample that is exactly in the middle between my measured chunks: $s(n=N+10000)$. I don't care about all the other samples in the gap.
What is the best and most efficient way to estimate my desired sample?
Note that $N$ would be ideally between 100 to 1000 and my sampling rate $f_s=100\mathrm{MHz}$.

One natural way to do this is to use the the total $2N$ samples (two chunks of $N$ samples each) to estimate the parameters. Amplitude, frequency, and phase estimation are all textbook equations that you can look up so the problem here is the fact that they are all mixed together. If you are assuming that you don't know how many sinusoids are mixed, then you also have to detect how many are present. The basic pseudo code is:
Detect how many sinusoids there are
For each sinusoid:
    Estimate amplitude
    Estimate frequency
    Estimate phase
end
Reconstruct signal at n = N + 10000

Detect how many sinusoids there are
You can do this by detecting how many peaks there are in the FFT of the signal. A sinusoid will present a peak at its frequency. At a peak, the slope of the curve will change from positive to negative, so the slope has a zero crossing at the peak. For a discrete signal, there are different ways to compute the slope/derivative but I would start with the simple before moving to more sophisticated methods. A simple way to compute the derivative is just the difference, $x_d[n]=x[n]-x[n-1]$, and a zero crossing happens if $x_d[n-1] > 0 > x_d[n]$. This will give you the number of zero crossings and where they occur.
Edit
The model that was presented in the question is a noise free model, and the detection described relies on that. You can still use this zero crossing method for a noisy signal but need to apply a threshold because the noise will cause many small peaks to occur. The threshold is needed to keep only the ones that have a significant enough amount of energy. Threshold should be picked to satisfy some probability of false alarm metric that you desire.
Estimate ...
As I mentioned above, estimating the amplitude, frequency, and phase of a sinusoid can be found through the FFT peak finding method (and we already found the locations of the peaks) so we are all set.
Reconstruct signal
Now that we have all the estimates, you can calculate your equation using your estimates.
$$ \hat{s}(n+10000) = \sum_i \hat{A}_i \text{sin}\bigg( \frac{2\pi \hat{f}_i (n+10000)}{f_s} + \hat{\theta}_i \bigg) $$


Where can I get an open source implementation of the Parks-McClellan optimal FIR filter design algorithm? Background: Often times I am doing some sort of signal processing task that requires a unique filter. Usually at this point I go to MATLAB and generate a new unique filter using $\tt firpm()$. The MATLAB firpm() function implements that Parks-McClellan algorithm. Now I have a filter and I put the filter into a hardcoded array. But here's the problem I now have a hardcoded filter that only works for one scenario.

The problem: I can now solve my signal processing problem du-jour... but only for a very SPECIFIC single sample rate or SPECIFIC scenario.

The goal: I want to be able to call $\tt firpm()$ from C code or some other language and make my signal processing code more generic. I can't find an open source implementation of firpm() !

Where can I get an open source implementation of the Parks-McClellan optimal FIR filter design algorithm (aka $\tt firpm()$ in MATLAB)?

P.S. I am aware that I can design filters differently using windowing or something else... feel free to mention those in the comments. But the point of this question is not to ask "what are other filter design techniques?" the point is to find an open source implementation of the VERY VERY useful $\tt firpm()$... or something similar.
P.P.S. One of the goals of this question is to learn how the Parks-McClellan algorithm works by looking at the code first and then I plan on reading some background theory.

There is an open-source implementation of Parks-McClellan (also known as the Remez exchange algorithm) in GNU Octave, a free-software implementation of a MATLAB-like environment. The function, called "remez", is contained in the "signal" package, which is hosted at Octave-Forge. If you download the package, you'll find "remez.cc", a C++ implementation of the algorithm.

One nice thing about Octave is that it is almost code-compatible with MATLAB, so you can easily port code over to use there if you like. It's a good way to peek under the hood at implementations of algorithms that are provided in MEX form in MATLAB.


Image noise reduction I have an image of a cross-section of soil where the main object of interest are the plant roots.

I would like to segment or extract the roots from the background noise of the soil itself and other obstructions such as spider webs. The goal is to form a clean image of the root network or end up with a result similar to what is shown here. In some cases, cracks in the soil itself are detected as interesting lines (structural artifacts) which is an unwanted result.
It seems I may have to apply a series of filters, does anyone have a suggestion what I could do to clean such image(s) or a pre-processing approach that would provide better results.
Current results

Maybe some of the following steps might help:

Desaturate the image using any method your prefer (either a color register, luminance or lightness). It is much easier to perform the next steps if the image is only 1 dimensional
possibly apply some non-linear function over the pixels gray value. This might increase the visual color difference between light pixels and dark pixels. Additionally, it will reduce the color difference between dark pixels and very light pixels. you could consider it as a smooth threshold function.
As the roots have a light color and the ground is overall dark, we should look for 2 aspects: a pixel with a light value and to which the neighboring pixels have a dark value. To enhance this second aspect, you could add a Laplace-edge detect filter to the image (so the resulting value will be the current value and the value of the laplace filter).
At this point, if properly done, the roots should already somewhat pop out the image. One could follow this up with some Dilate - Erode filters (you can stack those) to even out background color differences.

Hope this helped a bit.


Parallel equalizer, working for boosts but not cuts I am trying to replicate the behaviour of a parallel equalizer.

The approach is to use a parallel sum of bandpass filters like so:

$Response = 1 + bandpass1 * (gain1 - 1) + bandpass2 * (gain2 - 1) + ...$

Boosts and cuts are separately grouped and run in series.

I have successfully got the model working for two boosts in parallel. But running two cuts in parallel produces an unexpected result.

The figure below illustrates the problem. The thick blue line illustrates two boosts in parallel (thin orange and red line show the individual boosts that are being combined). The thick purple line illustrates two cuts in parallel (thin green and light-blue line show individual cuts that have been combined). As can be seen the individual cuts are a mirror image of the individual boosts but the parallel sum doesn't work for cuts.

Does anyone know what I'm missing so that that the thick purple line is a mirror image of the think blue line?

Here as my matlab code:

% boost 1
f1 = 176;
g1 = db2gain(12);
q1 = 1;

% boost 2
f2 = 672;
g2 = db2gain(12);
q2 = 1;

% cut 1
f3 = 176;
q3 = 1;
g3 = db2gain(-12);

% cut 2
f4 = 672;
g4 = db2gain(-12);
q4 = 1;

% bandpasses for boosts
[b1, a1] = bandpass(f1, proportionalQ(q1, g1));
[b2, a2] = bandpass(f2, proportionalQ(q2, g2));

% bandpasses for cuts
[b3, a3] = bandpass(f3, proportionalQ(q3, g3));
[b4, a4] = bandpass(f4, proportionalQ(q4, g4));

f = logspace(log10(20), log10(20000), 10000);

r1 = response(b1, a1, f); % bandpass response for boost 1
r2 = response(b2, a2, f); % bandpass response for boost 2
r3 = response(b3, a3, f); % bandpass response for cut 1
r4 = response(b4, a4, f); % bandpass response for cut 2

% run the analog bandpass filters in parallel as two boosts (OK)...
p1 = 1 + r1 .* (g1 - 1) + r2 .* (g2 - 1);

% run the analog bandpass filters in parallel as two cuts (NOT OK)...
p2 = 1 + r3 .* (g3 - 1) + r4 .* (g4 - 1);

figure;

semilogx(f, response2db(p1), 'LineWidth', 2); % parallel boosts (OK)

hold on;

semilogx(f, response2db(1 + r1 .* (g1 - 1))); % individual boost plot 1
semilogx(f, response2db(1 + r2 .* (g2 - 1))); % individual boost plot 2

semilogx(f, response2db(p2), 'LineWidth', 2); % parallel cuts (NOT OK)

semilogx(f, response2db(1 + r3 .* (g3 - 1))); % individual cut plot 1
semilogx(f, response2db(1 + r4 .* (g4 - 1))); % individual cut plot 2

xlabel('Frequency');
ylabel('Magnitude (dB)');

axis([20, 8000, -16, 16]);
grid;

hold off;

bandpass.m

function [b, a] = bandpass(f, q)
    b = [0, f / q, 0];
    a = [f * f, f / q, 1];
end

response2db.m

function [db] = response2db(r)
   db = 20 .* log10(abs(r));
end

db2gain.m

function [g] = db2gain(db)
    g = 10 .^ (db ./ 20);
end

proportionalQ.m

function [q] = proportionalQ(q, g)
    q = q .* sqrt(g);
end

response.m

function [r] = response(b, a, f)
    s = 1i .* f;

    r = (b(1) + s .* b(2) + s .* s .* b(3)) ...
        ./ (a(1) + s .* a(2) + s .* s .* a(3));
end

Cuts can simply be derived from boosts like so: p2 = 1 / p1

This is intuitive as this is how you would invert a bandshelf.

However, I'm not sure it leads to an elegant equation as inverting the sum of filters and recovering a rational filter expression might be messy, so there's potential a better way.


what does mean quasi-static channel When I read papers about channels, I usually see something Quasi-static Channel,

Does "Quasi-static Channel" mean "time-variant channel" or "time-invariant channel"?

Quasi-static channels can be said as "block-wise" time (in)variant. For example, your channel has no variations over time (not in delay domain) for some time-period, e.g., 1 ms, but the channel may change after 1 ms.


Convolving a discrete time sequence with a continuous-time filter? I am a newbie to DSP, and have what seems to me a fundamental question.

I have a discrete-time sequence $g[n]$ of length $N$. I am supposed to pass this signal through a root-raised cosine filter $h_1(t)$ with a cut-off frequency of $F= 1/T$, to obtain a continuous time signal $x(t)$:
$$
x(t) = \sum\limits_{n=0}^{N}g[nT]h_1((t-n)T)
$$

Is this a correct expression to obtain the continuous-time signal?
Then, $x(t)$ is convolved with a second continuous-time filter $h_2(t)$ to obtain $y(t)$:
$$
y(t) = x\star h_2 (t) = g\star h_1\star h_2(t)
$$

where $\star$ denotes convolution. I am not sure how to expand the above expression in terms of summation or integral, as $g[n]$ is discrete-time and the other two filters are continuous-time.
Finally, I sample $y(t)$ at rate $F = 1/T$ and pass it through $h_3[n]$ to obtain the signal $z[n]$. The filter $h_3[n]$ the discrete-time counterpart of $h_1(t)$, i.e. $h_3[n] = h_3[nT] = h_1[nT]$. Here, I have the same problem as 2. above. I am not sure how to expand the following in terms of summation or integrals, now that this is a mix-up of discrete-time signals $g[n]$, $h_3[n]$ and continuous time signals $h_1(t)$ and $h_2(t)$:
$$
z[n] = z(nT) = y(nT)\star h_3(nT) = x\star h_2\star h_3(nT) = g\star h_1\star h_2\star h_3(nT)
$$

Can I interchange the convolutions in the above as follows?
$$
g\star h_1\star h_2\star h_3(nT) = g\star h_1\star h_3\star h_2(nT)
$$

Any help or hints would be greatly appreciated!

-ryan

Conceptually, if used as an input to a continuous-time system, a discrete-time signal is represented as a weighted sum of Dirac delta impulses (as pointed out in a comment by MBaz):

$$g(t)=\sum_kg[k]\delta(t-kT)\tag{1}$$

If $g(t)$ in $(1)$ is the input to a continuous-time LTI system with impulse response $h(t)$, the output is given by

$$y(t)=\sum_kg[k]h(t-kT)\tag{2}$$

This expression is similar, but not identical, to the one you got. In your case, the impulse response $h(t)$ is simply the convolution of $h_1(t)$ and $h_2(t)$:

$$h(t)=(h_1\star h_2)(t)\tag{3}$$

Sampling $y(t)$ with a sampling interval $T$ means to replace the continuous time variable $t$ by $nT$ in $(2)$:

$$y[n]=y(nT)=\sum_kg[k]h((n-k)T)\tag{4}$$

If we define a discrete-time impulse response $h[n]=h(nT)$, we can write $(4)$ as a discrete-time convolution:

$$y[n]=(g\star h)[n]\tag{5}$$

Now we just have to filter $y[n]$ with the discrete-time impulse response $h_3[n]$ to obtain the final output signal:

$$\begin{align}z[n]&=(y\star h_3)[n]\\&=(g\star h\star h_3)[n]\end{align}\tag{6}$$

with

$$h[n]=h(nT)=(h_1\star h_2)(nT)=\int_{-\infty}^{\infty}h_1(\tau)h_2(nT-\tau)d\tau\tag{7}$$


how can I analize music as spotify/musixmath/shazam does? I mean, does exists some state of the art algorithm able to extract music features, classify them and identify similar music or find lyrics?

I know that there are many topics or steps involved on my question, like digital signal processing, feature extraction, machine learning, pattern recognition, etc. I'm a computer science student and I would like to know if there is a methodology to identify music as it does musixmatc. For example, this program is able to find the lyrics of a song from a very short audio sample, even with some noise. As I know, [1] is the unique related paper that is available.

1 An Industrial-Strength Audio Search Algorithm

The general idea for the algorithm is well presented in the linked paper. As is the case with applied machine learning and related fields, the secret sauce (which features to select, how to tweak your algorithms, etc.) is found by trial and error by people who know things about audio/audio signal processing (see for example, Julius O. Smith III's page at CCRMA @ Stanford) and by just trying a ton of things. Unfortunately, to become (and stay) the state of the art, the secret sauce has to be kept, well, secret.

That being said, the wikipedia article on acoustic fingerprinting is a decent place to look for some alternatives (and some links to some open source ones like Acoustid).


Relationship between Frequency response, magnitude spectrum and gain I have been reading quite a lot of DSP related materials online and offline (books etc.) lately. I come across the terms frequency response and gain. I have also noticed that both are expressed in dB.

I do understand the fact that expressing the magnitude spectrum in dB is more informative and intuitive and helps identify the presence of different frequencies. But as seen here (or here) and here, the plots basically have the same/identical parameters on X-axis as well as Y-axis i.e, its dB vs frequency. This makes me a bit confused.

What is the relationship between frequency response and gain. Are they the same thing? Are they used according to the context (in that case on what occasions it is apt to use either of them?)?

ok let me give you a practical (discrete-time) example.

Consider a system with input-output relation of the form
$$ y[n] = x[n] + 0.5  x[n-1] + 0.25 x[n-2] $$

This input-output relation satisfies linearity and time-invariance properties and hence defines an LTI (linear time invariant) system and therefore it has an impulse response $h[n]$ found to be :
$$h[n] = \delta[n] + 0.5 \delta[n-1] + 0.25 \delta[n-2] $$

Which signifies an FIR (finite impulse response) filter. Now since the system has an impulse response, then taking its Fourier transform will give us the system's so called complex valued Frequency Response function, denoted as $H(\omega)$ :

$$H(\omega) = \sum_n h[n] e^{-j \omega n} $$

For the above system, the frequency response can be shown to be:

$$ H(\omega) = 1 + 0.5 e^{-j \omega} + 0.25 e^{-j2\omega} $$

The magnitude of $H(\omega)$ signifies the (linear) gain that the LTI system applies to its input at the frequency $\omega$. For example, if $ |H(\omega)| = 0$ the system completely attenuates the input at that frequency, or if $H(\omega) = 1$ the signal passes unaltered and if $ |H(\omega)| = 10 $ then the system applies a linear gain of $10$ which corresponds to a dB scale of $20$ at that frequency.

Let's find the magnitude of the frequency response $H(\omega)$ :

$$
\begin{align}
H(\omega) &= 1 + 0.5 e^{-j \omega} + 0.25 e^{-j2\omega} \\
&= 1 + 0.5 (\cos(\omega)-j\sin(\omega) ) + 0.25 (\cos(2\omega)-j\sin(2\omega) ) \\
&= [ 1 + 0.5 \cos(\omega) + 0.25 \cos(2\omega)] -j [0.5 \sin(\omega) + 0.25 \sin(2\omega) ] \\
&= A(\omega) - j B(\omega) \\
|H(\omega)| &= \sqrt{ A^2(\omega) + B^2(\omega) } \\
|H(\omega)| &= \sqrt{ [ 1 + 0.5 \cos(\omega) + 0.25 \cos(2\omega)]^2 + [0.5 \sin(\omega) + 0.25 \sin(2\omega) ]^2 } \\
\end{align}
$$

with the help of a computer you can evaluate this frequency response magnitude at any frequency you want and compute the gain of that filter that it applies to any input frequency at that frequency. Below figure shows the Frequency Response magnitude of this system (evaluated and plotted using Matlab)

As you can see, the gain is about $1.8$ at frequencies close to zero radian per sample, and the gain is about unity at frequencies close to $0.5 \pi$ radians per sample. Most often but not always, this gain is defined in terms of a dB (deci Bell) scale (instead of a linear one) using the formula:

$$ \text{Gain_dB} = 20 \log_{10}( \text{Gain_linear}) = 20 \log_{10}( |H(\omega)|) $$

which is plotted in the figure below:

AS you can see, at the frequency about $0.5 \pi$ radians per sample, where the linear gain is about $1$, the corresponding deciBell gain is about $0$ dB.


How to calculate the energy in the dirac delta function signal? I'm having DSP for the first time, and after some classes I got confused about the following:

Suppose I have a signal which its fourier transform in a frequency band $[ \omega_1,\omega_2] $ is just a delta function with a peak somewhere in that band . What is the energy of the signal in that specific frequency band?

It seems that the answer depends on the integral

$$ \int_{\omega_1}^{\omega_2}|\delta (\omega-\omega^*)|^2.d\omega $$

with $\omega_1 < \omega^* < \omega_2$.

I'd guess this is 1 (maybe 0...), but I don't know how to compute it.

\begin{align}
&\int_{\omega_1}^{\omega_2}|\delta (\omega-\omega^*)|^2.d\omega \\
&= \int_{\omega_1}^{\omega_2}|\delta (\omega-\omega^*)|.|\delta (\omega-\omega^*)|.d\omega \\
&= \int_{\omega_1}^{\omega_2}\delta (\omega-\omega^*).\delta (\omega-\omega^*).d\omega \\
&= \int_{\omega_1}^{\omega_2}\delta (\omega^*-\omega^*).d\omega \\
&= \int_{\omega_1}^{\omega_2}\delta (0).d\omega \\
&= \delta (0).\Delta \omega
\end{align}
Which is undefined! See this for better clarification: https://en.wikipedia.org/wiki/Distribution_(mathematics)#Adding_and_multiplying_distributions

What you have to understand is that $\delta$ is not a function, it is a distribution. What you learn in introductory DSP courses is just an "engineering trick" so that things will work the way they should (you will understand what that means once you study sampling and analog to digital conversion).

EDIT: Sorry I misinterpretation the question. (Thanks Jazzmaniac for pointing that out)

As far as I know the $\delta$ "function" in the continuous domain is not in $L^2$, it isn't a square-integrable function.

Think about your question in the time domain. If you have 2 symmetric $\delta$'s in the frequency domain, you will have a INFINITLY lasting sinusoid. It's energy is infinite, i.e, you can't associate a value to it (undefined, remember that $\infty$ is not a number). If you have a single $\delta(f)$, centered at zero, you will have a complex exponential in the time domain, and it's energy is also undefined (unbounded since $|e^x|^2=1$.)


How to trace the edges on matlab I am doing my Bachelor's project on Road extraction. I have detected the line segments(edges) in the image using hough transform. I need to detect the curved lines in the road and link them(edge linking). help me with matlab coding.

One possible pre-processing apporach is to use imdilate(I,strel('disk',5)) to thicken the road lines a little bit.

If you are trying houghlines, probably one parameter that you need to tune is the 'FillGap' which will connect the very short segments (you can try to change this parameter yourself to observe the effect, along with the MinLength value):

lines = houghlines(BW,T,R,P,'FillGap',30,'MinLength',15);

However, the road curves are not straight lines, so I won't recommend you to use hough transform. You may try bwmorph(bw,'skel',N); after you thicken the binarized road image. N is the iteration number that you need to tune.


Independent component analysis for one observation of a Signal So I am a complete novice to ICA, so excuse my question if it is bad one, but I have the signal:

$$\sin(2\pi x) + \sin(4\pi x) + \textrm{Additive White Gaussian Noise}$$

I want to try to separate the two signals and apply a Fourier Transform on each one individually, to see if I can get a better spectral estimate than if they were together. I was thinking of using scikit's FastICA Algorithm. So my plan was to take the signal of length $3200$ samples, split it into $16$ windows each with length $200$ samples, apply a Hanning Function on each window, and do an ICA on 4 subsequent windows at a time.

Would this theoretically work, or does the ICA have to act on two different observations of a signal?
Is the ICA dependent on each observation of the signal having a different mixing matrix, or can it be the same?

You have single channel / observation along x. Splitting it along x would not make it multiple channels. I think you should not split it. You should rather approach it as underdetermined ICA problem, i.e., the number of sources are more than the number of observations.
For the point, theoretically, all the sources should be transformed via a system. So the observations, generated using different mixing matrices cannot be utilized jointly in a single ICA process.


Calculating the 3dB frequency of an nth order low pass filter and nth order high pass filter I am using the 6th order high pass filter in the C code. And the cut-off frequency selected for it is 150 Hz and the sampling frequency being used is 32000Hz. Now, I want to calculate the 3 dB frequency for this implementation (theoretically). I have read that the 3 dB frequency is the frequency at which the filter magnitude response is $1 / \sqrt{2}$ of the peak amplitude response. But I am not sure about how to calculate it? And also, I want to know how to calculate 3 dB frequency for an nth order low pass filter.
Note:
I am using the Chebyshev type II design in my code to implement a 6th order high pass filter.

In general there is no straightforward analytical solution. As you know, you need to solve
$$\left|H(e^{j\omega_c})\right|=\frac{1}{\sqrt{2}}\tag{1}$$
for $\omega_c$, where it is assumed that the maximum filter gain equals $1$.
For Butterworth filters, the specified cut-off frequency always equals the $3\textrm{ dB}$ frequency. This is not the case for other filter types.
In general, you have to either estimate $\omega_{c}$ from the magnitude plot, or you solve $(1)$ numerically.

EDIT: For Chebyshev filters it's relatively straightforward to compute an analytical expression for the $3\textrm{ dB}$ cut-off frequency.
As an example I'll use a Chebyshev II highpass filter, as mentioned in the question. The magnitude of the normalized frequency response of the analog prototype filter is given by
$$\left|H(\Omega)\right|=\frac{1}{\sqrt{1+\frac{1}{\epsilon^2T_n^2(\Omega)}}}\tag{2}$$
where $T_n(\Omega)$ is the $n^{th}$-order Chebyshev polynomial, and $\epsilon$ is a parameter determining the ripple size in the stop band. The normalized $3\textrm{ dB}$ cut-off frequency $\Omega_c$ can be determined from the equation $\left|H(\Omega_c)\right|=1/\sqrt{2}$, which is equivalent to
$$T_n(\Omega_c)=\frac{1}{\epsilon}\tag{3}$$
From the trigonometric definition of the Chebyshev polynomials
$$T_n(\Omega)=\cosh\big[n\; \textrm{arcosh}(\Omega)\big],\qquad\Omega>1\tag{4}$$
we can solve $(3)$ for $\Omega_c$:
$$\Omega_c=\cosh\left[\frac{1}{n}\; \textrm{arcosh}\left(\frac{1}{\epsilon}\right)\right]\tag{5}$$
Now we just need to transform $\Omega_c$ to the discrete-time domain according to the warping introduced by the bilinear transform:
$$\omega_{c}=2\arctan\left[\Omega_c\tan\left(\frac{\omega_s}{2}\right)\right]\tag{6}$$
where $\omega_s$ is the specified stop band edge of the discrete-time Chebyshev II high pass filter.
Note that we can also rewrite $(6)$ to determine the required value of $\omega_s$ given a specified $3\textrm{ dB}$ cut-off frequency $\omega_c$:
$$\omega_s=2\arctan\left[\frac{1}{\Omega_c}\tan\left(\frac{\omega_c}{2}\right)\right]\tag{7}$$
The following Matlab/Octave script illustrates the procedure of computing the $3\textrm{ dB}$ cut-off frequency of a Chebyshev II high pass filter:

% example filter specs of Chebyshev II high pass filter
n = 6;          % filter order
Rs = 60;        % stop band ripple in dB
ws = .5 * pi;   % stop band edge of Chebyshev II filter

% Cheby II: stopband ripple = ep / sqrt( 1 + ep^2 )
ep = 1 ./ sqrt( 10^( .1 * Rs ) - 1 );

% compute normalized 3dB cut-off frequency of analog prototype
Wc = cosh( acosh( 1 / ep ) / n );

% transform 3dB cut-off frequency to discrete-time domain
wc = 2 * atan( tan( .5 * ws ) * Wc );       % bilinear transform

I designed a discrete-time $6^{th}$-order Chebyshev II high pass filter according to above specifications. The resulting $3\textrm{ dB}$ cut-off frequency obtained from above script is $\omega_{c}=0.69\pi$. This is the correct value as shown in the plot below:


Applying estimated transfer function to a signal This is the simple code to find transfer function between sigout and sigin signals and then are the filter coefficients estimated by invfreqz function.

%% load simlated data
load CFDvsABAQUS_TC.mat

%% simulated data preprocessing (grid orientation B)
time = TCb(:,1);
sigin = TCb(:,2);           % inlet temperature (TC signal by Abaqus)
sigout = TCb(:,3);          % outlet temperature (by cfd FENIX)
Fs = 1/mean(diff(time));    % sampling frequency

%% Transfer Function Estimate outlet->inlet
[toi,foi]=tfestimate(sigout,sigin,2^13,[],[],Fs);

%% filter design by known transfer function
w = linspace(0,pi,length(fio)); % angular frequency f=[0,fs/2] -> w=[0,pi]
m = 500; % b coeffs
n = 500; % a coeffs
[boi,aoi] = invfreqz(toi,w,n,m);

%% final Digital filtering
esigin = filter(boi,aoi,sigout);

The input data are available at this dropbox link.

The final filter function produce completely wrong results (the result should be similar to original sigin signal). I have no idea what is wrong.

Any help will be very useful for me.

Your values for m and n are way too big. Matlab starts getting numerical problems for those values above 30, sometimes.

Try reducing them to more reasonable numbers (e.g. 10) and slowly increase them until invfreqz generates bogus outputs.


Nyquist Theorem adding two same frequency near to Nyquist Frequency with phase shift This is my first question on this platform. Sorry if I made mistakes.

What happens if we add two or more same frequency signals near to Nyquist Frequency with phase shift and sample them?

For example, assume that we have a signal of 18000Hz. And add another signal with same frequency, but having 180 degree phase with our first signal. Then, sample this total signal at 36000Hz. What happens if we reconstruct this signal? Does this situation result in aliasing?

Note: I'am trying to figure out how exactly sound records are sampled when musicians playing their instruments. For example, assume that two keyboard players playing a note at 20000Hz which will be sampled at 44100Hz. Since they are human-being they are not supposed to hold the rhythm like in electronic music. Therefore, there will be a phase shift in their playing and wouldn't this cause aliasing in reconstructed signal?

Thanks for the answers.

When you add 2 or more sinusoids at the same frequency $f_o$ but with different phase shifts you get a sinusoid at same frequency but an additional attenuation term. Mathematically, you will have following:
$$cos(2\pi f_ot + \phi)+cos(2\pi f_ot + \theta) = 2cos(\frac{\phi - \theta}{2}).cos(\frac{2\pi f_ot + \phi + 2\pi f_ot + \theta}{2})$$$$= 2cos(\frac{\phi - \theta}{2}).cos(2\pi f_ot + \frac{\phi + \theta}{2}),$$You see basically, you have an attenuation term $2cos(\frac{\phi-\theta}{2})$ which will change the amplitude of the sum of the sinusoids, but the frequency of the sinusoid remains at $f_o$. Hence, if you were sampling it with $f_s>2f_o$, then even after adding them you will not have any aliasing.

Only thing of concern is that because of attenuation term driving the amplitude of the sinusoidal, your signal might become 0 now, as in you example of $2^{nd}$ sinusoidal at the same frequency $f_o$ but having phase difference($\phi - \theta$) of $\pi$ compared to $1^{st}$ signal. That means $2cos(\frac{\pi}{2})$ which is equal to 0, and hence you have actually nothing to sample, both signals are completely out of phase and cancelling each other.


How to perform a Rubberband-Correction on spectroscopic data? I basicaly understood how the Rubberband/Baseline-correction works.

The given spectrum is divided into (N) ranges.
The lowest points in every range are determined.
The initial baseline is built out of those points.
Now all the points on the spectrum are drawn down by the difference between the lowest point in the current range and the lowest point on the baseline.

There are some nuances, though, that I do not know how to handle. E.g., what if one of the points is exactly on the border between two ranges, etc.

Plus, I have to be able to prove that the algorithm that I am writing is a solid one and can be referenced by other works or scientific papers.

If anyone could give me some reference I would be very pleased.

This can be easily done in R or Python. There are well-tested functions available, so you don't have to worry about any boundaries or nuances. Moreover, both are free and popular among scientists.

Solution for R

There is a special package to handle spectral data, called hyperSpec. The rubberband baseline correction is already implemented there (function spc.rubberband). All details are highlighted in the documentation. The usage is pretty simple:

require(hyperSpec)
spc <- read.spe("paracetamol.SPE")
baseline <- spc.rubberband(spc)

corrected <- spc - baseline

Solution for Python

There is no (to the best of my knowledge) out-of-the box solution for python, but you can use scipy.spatial.ConvexHull function to find indices of all points that form a convex hull around your spectrum. Suppose that the spectrum is contained in x and y arrays:

import numpy as np
from scipy.spatial import ConvexHull

def rubberband(x, y):
    # Find the convex hull
    v = ConvexHull(np.array(zip(x, y))).vertices

Array v contains indices of the vertex points, arranged in the CCW direction, e.g. [892, 125, 93, 0, 4, 89, 701, 1023]. We have to extract part where v is ascending, e.g. 0–1023.

    # Rotate convex hull vertices until they start from the lowest one
    v = np.roll(v, -v.argmin())
    # Leave only the ascending part
    v = v[:v.argmax()]

    # Create baseline using linear interpolation between vertices
    return np.interp(x, x[v], y[v])

Now the baseline is corrected like this:

y = y - rubberband(x, y)


what is the theory behind Butter function? I want to understand why this code used tan and cos ?

function [ b, a ] = butterTwoBp( dt, fl, fu )
q=pi*dt*(fu-fl);
r=pi*dt*(fu+fl);
N = (tan(q)^2) + sqrt(2)*tan(q) + 1;
M = (tan(q)^2) / N; %M after N because it depends on N
O = -cos(r) * (2*sqrt(2)*tan(q) + 4) / ((cos(q))*N);
P = (-2*(tan(q)^2) + ((  (2*cos(r))   /  (cos(q))   )^2) + 2 )  /   N;
Q = cos(r)*(2*sqrt(2)*tan(q) - 4)/(cos(q)*N);
R = (   (tan(q)^2) - sqrt(2)*tan(q) + 1   )  /  N;

b=[M 0 -2*M 0 M];
a=[1 O P Q R];

I got this code from https://stackoverflow.com/questions/17340198/what-is-the-command-for-butterworth-bandpass-filter?answertab=active#tab-top

The trigonometric functions enter the calculation because of the use the bilinear transform for transforming an analog filter to a discrete-time filter. The bilinear transform warps the frequencies of the analog filter, that's why we have to pre-warp the frequencies of the analog filter in order to obtain the desired cut-off frequencies of the discrete-time filter.

The frequency warping of the bilinear transform is described by the following equation:

$$\Omega=\tan\left(\frac{\omega}{2}\right)\tag{1}$$

where $\Omega$ is the frequency (in radians) in the analog domain, and $\omega$ is the normalized frequency (in radians) in the discrete-time domain. Eq. $(1)$ is used to pre-warp the specified frequencies to obtain the edge frequencies for designing the analog prototype filter. After applying the bilinear transform to that analog prototype, the resulting discrete-time filter will have the correct edge frequencies.

The design procedure is as follows:

Pre-warp the given edge frequencies according to $(1)$.
Design an analog band pass filter with the pre-warped edge frequencies by applying a frequency transformation to a low pass prototype filter.
Apply the bilinear transform to the analog band pass filter to obtain the desired discrete-time filter.


Communication between two USRP : mismatched spectra between transmission and reception I am trying to make two USRP communicate between each other : for the transmission, I'm using an HackRF One USRP and for the reception a NI 2921 USRP.
As a first test, I'm trying to send different basic signals like sine waves, square waves, triangle waves,... with the HackRF One and receive it with the NI 2921 USRP at a carrier frequency of 2.43GHz.

Here is the block diagram for transmission, the time signal and FFT plot :

I'm sending the signal at 2.43GHz with a sample rate of 200kHz. Here the signal sent is a composition of a real component which is a square wave at 15kHz and an imaginary component which is a cosine wave at 10kHz.
Here is the block diagram for reception, the time signal and FFT plot (once when I am not sending any signal and once when I am transmitting the composition of signals described above, such that we can see that the peak detected is well due to my emission) :

I'm sampling at 200kHz and my center frequency is at 2.43GHz.

Whatever the signal I'm sending (a DC signal, a sine wave, a triangle wave,...), I always get the same peak (and only one) like the one on the picture. Normally, I should receive the same FFT as the one showed at emission on the picture.
I notice also that this peak is located at a weird position (2.43GHz + 46kHz) compared to the center frequency (2.43GHz) with no link with the frequency of the signals sent. Moreover whatever the frequency of the signal emitted, the peak will always be located at this same position as well.
I seems like they manage to communicate between each other given that there is this peak appearing when I am transmitting the signals but I do not understand the logic behind the spectrum received..

Can somebody enlighten me?

EDIT further to @MBaz answer:

I tried both osmocom and UHD at the reception of the signal to see if there was a difference but there is not. I will use UHD again.
I have just used GQRX as a spectrum analyzer to do some tests. Sending a sine wave at 10kHz gave me peaks at frequencies shifted of 50kHz on the FFT (if I send a signal at 2,429720GHz, the received signal showed a peak at 2,429770GHz on the FFT with GQRX). I am then sure the received peaks are due to the emission of signals with the HackRF One.
Changing the sine frequency from 10kHz to 20kHz does not change the location of the peak.
Trying to send two sine waves (10kHz and 20kHz) showed me one and only one peak at the same frequency location as before.
Changing the type of signal sent (triangle, square wave, sine wave) does not change the FFT of the signal received.
All the previous tests have showed the same results on gqrx and GNURadio Companion.
As you can see on the plot of the evolution of the signal in function of the time, this is a quite unusual low amplitude signal received.
No error is shown on the console of the NI 2921 receiver. On the emitting signal console, I have this output :

Warning: failed to XInitThreads()
linux; GNU C++ version 7.3.0; Boost_106501; UHD_003.010.003.000-0-unknown
gr-osmosdr 0.1.4 (0.1.4) gnuradio 3.7.11
built-in sink types: uhd hackrf bladerf soapy redpitaya freesrp file
Cannot connect to server socket err = No such file or directory
Cannot connect to server request channel
jack server is not running or cannot be started
JackShmReadWritePtr::~JackShmReadWritePtr - Init not done for -1, skipping unlock
JackShmReadWritePtr::~JackShmReadWritePtr - Init not done for -1, skipping unlock
Using HackRF One with firmware 2018.01.1

A few (hopefully useful) comments and ideas:

The HackRF One is not a USRP.
If you're receiving with an NI 2921 USRP, you should be using UHD to interface with it, not osmocom.
Use a standalone spectrum analyzer, first to double-check that nobody else is transmiting at that frequency, and then to study your transmitted signal.
Transmit a single sine wave and verify that both the spectrum analyzer and the USRP can see it. Verify that you're transmitting it with enough power.
If that works, try with two sine waves separated by a few kHz.
It is normal that the received signals are at a slightly wrong frequency. Each device has its own idea of what 1 Hz is; this is why receivers need to run a carrier synchronization algorithm.
Your devices may not support a sampling frequency of exactly 200 kHz and may be using the closest available frequency. Check any warnings on the console. This may result in impaired performance.
In general, pay attention to any output on the console: you may be overrunning or underrunning your flowgraphs, or something else may be wrong.


How to get spectrum of an image? I have an image :

img = imread('pic.jpg');

Now, How can I get spectrum of that image?

reference : Digital Image processing - Rafael Gonzalez

Steps :
The whole of filtering process can be summarized in the following points :
1> Given an input image f(x,y) of size M X N, obtain the padding parameters P and Q as P =2M and Q = 2N

The problem with using DFT is that both the input and output images are periodic. This causes interference between adjacent periods if the periods are close with respect to the duration of the non-zero parts of the functions. So we padd  the image by using a padding function P having zeros equal to twice the image length.

2> Form a padded image f’(x,y) of size P X Q by appending the necessary numbers of zeros to f(x,y).

3> Multiply f’(x,y) by (-1)^x+y to center its transform

4> Compute the DFT, F(u,v) of the image.

f = imread('pic.jpg');
figure(1);
imshow(f);
title('original image');

PQ = paddedsize(size(f));
[U,V] = dftuv(PQ(1),PQ(2));
 F = fftshift(fft2(f,PQ(1),PQ(2)));

 S =  abs(F);
 ag = angle(F);

[l,p] = size(S);
d = max(max(S));
c = min(min(S));

for l = 1:1:l
    for p = 1:1:p
        S1(l,p) = (S(l,p)-c)*(255/(d-c));
    end
end

figure(2);
subplot(2,2,1);
imshow(uint8(S1));
title('Magnitude response without log transform');

% since dc dominates the values of spectrume, the dynamic range of other intensities in displayed image is compression. to bring out those details we perform log transformation.

[m,n] = size(ag);
d = max(max(ag));
c = min(min(ag));

for l = 1:1:l
    for p = 1:1:p
        ag1(l,p) = (ag(l,p)-c)*(255/(d-c));
    end
end

subplot(2,2,2);
imshow(uint8(ag1));
title('Phase response without log transform');


Non Linear Remapping of the Pixel Values Suppose you take a photo of a beautiful scene of mountains and a lake. The representation devotes half of the scale between black and white to the bright surface of the lake, and the other half to the subtle shadings of the mountains. Is it possible to change the mapping so that only ¼ of the scale is devoted to the lake's colors, and ¾ to the mountains shades? If so, how would you do it?

Yes, using a non-linear mapping function on the intensity. Let us for the moment suppose you have a grayscale image (because of the words "black and white", "shadings"), with luminance $0\leq p_{i,j} <1$.

The piecewise linear mapping: $f( p_{i,j}) = \frac{p_{i,j}}{2}$ if $p_{i,j}<\frac{1}{2}$, and  $f( p_{i,j}) = \frac{3}{2} (p_{i,j}-\frac{1}{2})+\frac{1}{4}$ if $p_{i,j}\geq\frac{1}{2}$ can do the job.

This is a simple example of "Non-linear Contrast Enhancement". If that works for you, you can try with smoother functions.


Z transform getting different answer for transform of rational function Why am i getting different answers and which one is correct??

In this also

Your first line is wrong. The partial fraction expansion of $X(z)$ should be :

$$ X(z) = \frac{ z (2z-\frac{5}{6})}{ (z-\frac{1}{2})(z-\frac{1}{3})} = 2 + \frac{ \frac{1}{2} }{z -  \frac{1}{2}} + \frac{ \frac{1}{3}}{z -  \frac{1}{3}} $$

Then they should yield the same inverse (causal) transform as:

$$ x[n] = 2\delta[n] + (1/2)^{n} u[n-1] + (1/3)^{n} u[n-1] $$

which is identical to :
$$ x[n] = (1/2)^n u[n] + (1/3)^n u[n] $$

as obtained by the second method.


Image 2D Real Cepstrum with DFT, Is `ifftshift` Needed? This is my testing image,it is taken from this paper:

I tried to transform it into its real cepstrum domain with this simple MATLAB code:

cepstrum_img=ifft2(log(1+abs(fft2(img(:,:,1)))));
imshow((abs(cepstrum_img)))

But my result is very weird: The maximum is not at the center of image, they are at the 4 corners. And here is the zoomed-in picture with the  maximum at the left top corner:

But the result in this paper is like this:

The maximum is actually at the center!!

As I didn't do any fftshift ,  I don't think I need to do any ifftshift.

This makes me so confused. Is there something wrong with my code?
or my understanding about 2D real cepstrum with DFT?

UPDATE 1: Why do I think that fftshift or ifftshift is not needed here.

For a length-$N$ sequence $x[n]$, you do an $N$-point DFT. You will get another length-$N$ sequence $X[n]$. The first element of the resulting sequence is corresponded to the zero frequency. Then we do fftshift to the resulting length-$N$ sequence, and we get a third length-$N$ sequence $X'[n]$, because we want to make the center (may not be at the center if $N$ is even) of the element in $X'[n]$ corresponded to the zero frequency.

The question is how can you get the original  $x[n]$ with IDFT?

$\textrm{IDFT}(X[n])$ is right. $\textrm{IDFT}(X'[n])$ is wrong.

Check this MATLAB code:

% the original sequence
x=1:8
% this will give you the original sequence
ifft(fft(x))
% but this won't
ifft(fftshift(fft(x)))
% anyway, the amplitude is the same.

The case would be the same for images and fft2 & ifft2. i.e.this will get you the original img:

 ifft2(fft2(img))

And if the cepstrum is defined as

$$C=\mathcal F^{-1}\left\{\ln\left(\lvert \mathcal F\left\{I(x,y)\right\}\rvert\right)\right\}$$

Where $I$ = img. This is something like transforming from Fourier domain back to the spatial domain.

UPDATE 2: I found where I am wrong..

When the papers are saying

  The Cepstrum (Calculated with FT) has the maximum at the center

I misunderstood it as

  The Cepstrum calculated with DFT & IDFT has the maximum at the
  center

If the Cepstrum above is really calculated with Fourier transform, the Cepstrum will be with the maximum at the origin. That means if I calculate this with DFT & IDFT,the maximum should be at the corners. The same as FT vs DFT, I have assumed the Cepstrum calculated with DFT should have its maximums at the center. That's where I am wrong.

If you have image in the spatial domain, in order to calculate its DFT transform you should use fft.
Once you have its DFT in order to get back to the spatial domain use the function ifft.

Either way, the DFT transform gives you the data in the [0, 2pi] axis.
Use fftshift to move it into the [-pi, pi] domain as seen in the paper.

Good Luck.

MATLAB Code

This code worked for me:

mInputImage = imread('3.png');
mInputImage = double(mInputImage) / 255;
mInputImage = mean(mInputImage, 3);

mImageDft = (fft2(mInputImage));

mImageDft = log(1 + abs(mImageDft));

mImageCepstrum = abs(fftshift(ifft2(mImageDft)));

mImageCepstrum = mImageCepstrum / max(mImageCepstrum(:));

figure();
imshow(mImageCepstrum);

Enjoy...


Relation Between Compressed Sensing and Incoherence I am working on CS. I know that as coherence decreases (or incoherence increases) compressed sensing results improve. I know how to calculate the coherence, but I can't find an exact equation between incoherence and compressed sensing. I would like to know the equation relating coherence and CS with lower and higher limits.

Thanks in advance.

The coherence can be used to bound the number of measurements $m$ required for successful reconstruction as seen for example in equation (6) in Candès & Wakin 2008:
$$m ≥ C \cdot \mu^2(\Phi, \Psi)\cdot S  \cdot \log n$$
where $C$ is some positive constant, $S$ is the sparsity, and $n$ is the size of the sparse signal.

One problem, however, is that this is not a very sharp bound. Donoho & Tanner 2010, for example, describe some brief history of this in Section X of the linked paper. I have not followed the development on this aspect of compressed sensing for a long time, so I do not know if there are newer or better explanations. So, in my opinion, the best you can use the coherence for in practice is as a rule of thumb: the less coherence, the better.


Good way for segmenting this signal I want to segment the following signal (units are x: [px], y: [a.u.]) and need a good way to do so. I want to achieve the following:

I am only interested in the two wide "plateaus" with maximum value (e.g. (1) from ~400px to ~800px and (2) from ~1300px to ~1800px).
Ideally would like the solution to be parameterless.

In the figure the blue line is the original signal and the orange line is the signal filtered by savgol_filter(signal, 51, 1).

I have tried Otsu's method, but it is not optimal, because it is intended for separating a foreground from the background. This is not the case here, which is why I had to do the smoothing (which introduces a parameter(!)) to reduce the noise in the side-peaks, so that they do not "leak" into the foreground, when using Otsu.

I have considered using the derivative of the signal, which should be very strong the flanks of the plateaus (see second figure), but here again the noise in the side-plateaus is causing an issue.

So what would be a good way to segment this signal?

Signal to segment:

signal_filtered = savgol_filter(signal, 51, 1)
plt.plot(signal)
plt.plot(signal_filtered)
plt.show()

Differential of signal and smoothed signal:

plt.plot(np.abs(np.diff(signal))), plt.plot(np.abs(np.diff(savgol_filter(signal, 51, 1)))), plt.show()

EDIT:

Signal with single, non-flat plateau:

I solved this by still using Otsu to find the threshold value. I then just filter the the regions with value == 1 in the thresholded signal (i.e. the region_mask; green line in plot below) by their width.
This is the result, where the rectangular peak on the right will be filtered (see below for the code):

plt.plot(signal_orig)
plt.plot(signal)
plt.plot(region_mask * np.max(signal))
plt.show()

Code:

signal = savgol_filter(signal_orig, 51, 1)  # window size 51, polynomial order 1
threshold = threshold_otsu(signal)

region_mask = signal > threshold
region_list = get_regions_from_mask(region_mask)
region_list = filter_date_regions_by_width(region_list, minimum_region_width=200)

def filter_date_regions_by_width(region_list, minimum_region_width):
    """
    Filters removes data regions from region_list with width < minimum_region_width.

    :param region_list: list of DataRegion objects, that we want to filter.
    :param minimum_region_width: minimum width of a region, so that it is not filtered out.
    :return:
    """
    filtered_list = []
    for region in region_list:
    if region.width >= minimum_region_width:
        filtered_list.append(region)
    return filtered_list

def get_regions_from_mask(region_mask):
    """
    Get list of the DataRegion objects that were found inside region_mask.

    :param region_mask: the np.array of 0 and 1 in which we want find regions of value == 1.
    :return: region_list: a list of DataRegion objects.
    """
    region_list = []  # list that will hold all regions that were found
    region = DataRegion()  # object to hold the start and end values of the region
    inside_region = False  # indicates, if we are currently inside a mask_region with value=1
    for index, value in enumerate(region_mask):
    if value == 1 and not inside_region:
        # entered a region
        inside_region = True
        region.start = index
    if value == 0 and inside_region:
        # left a region
        inside_region = False
        region.end = index
        region.width = region.end - region.start
        region_list.append(region)
        region = DataRegion()
    return region_list

Definition of the DataRegion class:

from dataclasses import dataclass

@dataclass
class DataRegion:
    start: int = None
    end: int = None
    width: int = None


When listening in to an AM signals of various frequencies, how do we exactly tune in? Alright, so I read on AM and Double Sideband AM, but I don't the more fundamental idea - the frequencies.

Frequency is just: $ \frac{1}{T} $, where T is period

So difference between 5 Hz and 30 Hz is that period of 1st signal is $\frac{1}{5}$ = 0.2 seconds, and of 2nd signal = $\frac{1}{30} = 0.03 seconds$

and I know in AM you need to demodulate by multiplying by cos(wt) to bring message to lower frequency and then low pass filter to only have message signal etc

but for simplicity sake, let's assume two signals are transmitted, one is at 5Hz the other one is at 30Hz

How do you "listen in"? Because the frequency spectrum, where frequency is x-axis is a bit misleading!

It shows that you can magically shift to the right and just listen in to the needed signal! That's not the case! Frequency is just how fast a signal repeats.
Amplitude is changing according to sine function.

So if you have two signals, one at 5Hz and second one at 30Hz, then they both "are in the air", you can't just magically listen in every $\frac{1}{30Hz} = 0.03 seconds$ and only get signal at 30Hz, you'll also be getting the signal at 5Hz, they'll be overlapping because frequency is just how fast signal is changing.

So how do you listen in? Do you just listen in for everything?

Then you have a signal that's changing every "0.2 seconds (5Hz)" and another one every "0.03 seconds", they'll overlap of course. The only way they don't is if there's only one signal is in the air. I don't see how different frequencies allow to fine tune to one particular signal and ignore others.

Besides, even if there's just one signal in the air, say it's the 5Hz signal, how do you "listen in"? Do you just "listen in" continuously, or a receiver just turns on and shuts off every $\frac{1}{5Hz} = 0.2 sec$????

I watched videos, and they just say "receiver" receives, but I don't get how it receives...

Ok, I think here is a good point to introduce equivalent baseband, because you implicitly are already using it!

So, what does your cosine 5Hz audio signal look like, if you were to set $f_c=0$, in spectrum?

Exactly, one dirac at +5 Hz, and one at -5 Hz! Hence, when you mix that up to $f_c\ne 0$, you get symmetric sidebands. That works for any real signals (ie. for any audio signal) – and even when you add them. That's where your two sidebands come from – positive and (hermitian) symmetrical negative components of any real-valued signal.

So, in baseband, your sum of 5 and 30 Hz cosines have four spectral components – at -30, -5, +5 and +30 Hz, and mixed up to the carrier frequency you get the same discrete spectral components, but added to $\pm f_c$. (your figure also shows a component at $f_c$, but that's only there if you got a DC component in your baseband signal, and/or you're not suppressing the carrier)

Now, there's very different receiver architectures for AM-modulated audio. None of them "just looks every 1/(audio period)". "Looking every so often" presumes you're building something digital. Most AM receivers (most AM being pretty obsolete by now) are not digital. The simplest detector methods really are just rectifying diodes and a low pass filter to get rid of the RF content – I'll leave googling for "AM diode detector" up to you; there's a plethora of good articles on that out there. This is all continuous-time, so there's no "looking that often" there – it all happens by electronically processing the continuous signal.

Now, assuming we're really aiming for digital here:

That 1/(audio period) wouldn't even make sense for audio – it breaks Nyquist; you need at least twice the audio bandwidth as sampling rate to be able to reconstruct the signal.

What one can build is simply a mixer with $f_c$, which will first (continuous-time!) multiply with a harmonic of that frequency (effectively, a complex sinusoid), and then sample (that's the act of looking only every so often) that, and then you get a digital signal that's basically nothing but your original audio sum.

Note that you've been asking about "how to deal with the fact that you add two sines of different frequencies", but that question is totally independent from the AM aspect of that: you just can. Under the Fourier Transform, i.e. in the spectrum, these are orthogonal, i.e. you can perfectly separate them, both in a computer, in a circuit, or with your ear. That's why you're hopefully able to listen to sounds that aren't made of a single tone. In fact, in reality, single-tone sounds are extremely rare.


Determine stability of feedback system from open loop transfer function and Nyquist stability criterion gives different results I'm confused due to the fact that the Nyquist stability criterion and looking at the transfer function doesn't give the same results whether a feedback system is stable or not. When I have the system of this block diagram the closed loop transfer function is:
$$
T_{CL}(s) = \frac{G}{1+GH}
$$
Due to the Nyquist stability criterion this closed loop transfer function will be marginally stable when GH = -1 + 0j and unstable when $GH = < -1 + 0j$ (if the open loop transfer function is stable, in this case there aren't any poles in the right half plane thus the nyquist plot should not encircle the point -1, so this plot should stay right from this -1 point).
However when I look at the closed loop transfer function, I would say that this system is unstable for $GH = -1$. In this case the transfer function becomes infinity so a bounded input will result in a unboundend (=infinity) output.
In my train of thought the point $GH = -2$ would again be stable since $T_{CL}$ will be finite again, however conform the Nyquist stability criterion this point wil still be unstable?
I know that Nyquist is correct but what's the problem with my way of thinking

However when I look at the closed loop transfer function, I would say that this system is unstable for =−1. In this case the transfer function becomes infinity so a bounded input will result in a unbounded (=infinity) output.

This depends on your definition of stability.  $GH = -1$ is called marginally stable because depending on how you look at it, it could be stable, or it could be unstable.
In the Lyapunov sense that @Petrus1904 mentions, it is stable.  But if you take bounded-input bounded-output (BIBO) stability to mean that the input can go on infinitely long but within bounds, then the output can, indeed, go to infinity.  So in the BIBO sense it's unstable.
Hence, "marginally stable".

In my train of thought the point =−2 would again be stable since  will be finite again, however conform the Nyquist stability criterion this point will still be unstable?

An unstable system can have a transfer function that never goes infinite for $s \in j \omega$.  All it requires is that there are no poles on the stability boundary.
I apologize for not taking the time to do so, but I know that I could dummy up a system that has $GH = -2$ and was still stable, because I've designed such systems.  All you need to do is wrap a double integrator with a PID controller: let $G = 1/s^2$ and $H = k_i / s + k_p + \frac{s}{\tau_d s + 1} k_d$ and tune for stability.  The resulting system will have a spot in its open-loop Bode plot which has a phase shift of 180$^\circ$ and a gain greater than one -- to get that gain equal to two you'll just need to jigger the numbers around.
Addendum:  Use the PID transfer function I give above with $k_i = 0.002$, $k_p = 0.02$, $k_d = 0.2$ and $\tau_d = 0.01$ and you should get a system that's stable, with $GH \simeq -2$ at somewhere around $\omega = 0.1 \mathrm{\frac{rad}{sec}}$.
Assuming I got my math right...


What is correlation kernel and compare with gaussian kernel I read a paper that said about correlation kernel that defined:
$$W(x-y)=(α/1+d(|y − x|))$$
where $α =  (\int(1+d(y − x)dy)^{-1}$, $(d(|y − x|))$  is spatial Euclidean distance from the central pixel.
My question is what is correlation kernel and benifit of that kernel comparison with gaussian kernel? Thank you so much.

Referencing the paper, the convolution kernel can be written as:

$d(|y-x|)$ is the Euclidean distance between the centre pixel $y$ and a location in the kernel $x$.  $\alpha$ is used to ensure that the entire area under the kernel is 1.  However, you did not specify how big this kernel is.  As such, we will specify the rows and columns of this kernel to be M and N respectively.  Let's also assume that the size of the kernel for each dimension is odd.  The reason why is because the shape of the kernel will be an even square and makes implementation easier.  As such, here are the steps that I would perform to do this:

Define a grid of X and Y co-ordinates, and ensure that the centre pixel is at 0.
Compute each term in the convolution kernel without the $\alpha$ term.
Sum up all of the terms in the kernel, then divide every value in this kernel by this term so that the entire area of the kernel is 1.

Let's do this step by step:

Step #1

We can do this by using meshgrid.  meshgrid (in this case) creates a 2D grid of (X,Y) co-ordinates.  X defines the horizontal co-ordinate for each location in X, while Y defines this vertically. By calling meshgrid(1:m, 1:n), I am creating a n x m grid for both X and Y, where each row of X progresses from 1 to m, while each column of Y progresses from 1 to n.  Therefore, these will both be n x m matrices.  Calling the above with m = 4 and n = 4 computes:

m = 4;
n = 4;
[X,Y] = meshgrid(1:m, 1:n)

X =

 1     2     3     4
 1     2     3     4
 1     2     3     4
 1     2     3     4

Y =

 1     1     1     1
 2     2     2     2
 3     3     3     3
 4     4     4     4

As such, we simply have to modify this, but we ensure that the middle is at (0,0), and also ensure that the size of X and Y are odd.  Let's say that M = 5 and N = 5.  We can then define our X and Y co-ordinates like so:

M = 5;
N = 5;
[X,Y] = meshgrid(-floor(N/2):floor(N/2), -floor(M/2):floor(M/2))

X =

-2    -1     0     1     2
-2    -1     0     1     2
-2    -1     0     1     2
-2    -1     0     1     2
-2    -1     0     1     2

Y =

-2    -2    -2    -2    -2
-1    -1    -1    -1    -1
 0     0     0     0     0
 1     1     1     1     1
 2     2     2     2     2

As you can see here, the centre pixel for both X and Y are defined as (0,0).  Everywhere else has its (X,Y) co-ordinates defined with respect to the centre.

Step #2

We simply have to compute the Euclidean distance between the centre pixel to every point in the kernel.  This can be done by:

dis = sqrt(X.^2 + Y.^2)

dis =

2.8284    2.2361    2.0000    2.2361    2.8284
2.2361    1.4142    1.0000    1.4142    2.2361
2.0000    1.0000         0    1.0000    2.0000
2.2361    1.4142    1.0000    1.4142    2.2361
2.8284    2.2361    2.0000    2.2361    2.8284

Doing some quick calculation checks, you can see that this agrees with our understanding of Euclidean distance.  Moving to the left by 1 from the centre is a distance of 1.  Moving to the left by 1 then up by 1 gives us a Euclidean distance of $\sqrt{1^2 + 1^2} = \sqrt{2} = 1.4142$.  Doing similar checks with each element will demonstrate that this is indeed a Euclidean distance field from the centre pixel.  After we do this, let's compute the kernel terms without the $\alpha$ term.

kern = 1 ./ (1 + dis)

kern =

0.2612    0.3090    0.3333    0.3090    0.2612
0.3090    0.4142    0.5000    0.4142    0.3090
0.3333    0.5000    1.0000    0.5000    0.3333
0.3090    0.4142    0.5000    0.4142    0.3090
0.2612    0.3090    0.3333    0.3090    0.2612

Step #3

The last step we need is to normalize the mask so that the total sum of the kernel is 1.  This can simply be done by:

kernFinal = kern / sum(kern(:))

kernFinal =

0.0275    0.0325    0.0351    0.0325    0.0275
0.0325    0.0436    0.0526    0.0436    0.0325
0.0351    0.0526    0.1052    0.0526    0.0351
0.0325    0.0436    0.0526    0.0436    0.0325
0.0275    0.0325    0.0351    0.0325    0.0275

This should finally give you the correlation kernel that you are seeking.  You can now use this in convolution (i.e. using imfilter or conv2).

Hopefully I have answered your question adequately.  Good luck!


How to undo the quantization step in DCT? How can I do point by point multiplication of the quantized DCT coefficients with the quantization matrix before computing the inverse DCT for each block? I want to modify the reconstruction algorithm to approximately undo the quantization step shown here:

for i = 1:Blocks
    for j=1:Blocks
        irange = (i-1)*8+1:(i-1)*8+8;
        jrange = (j-1)*8+1:(j-1)*8+8;
        quantblock = round(y(irange,jrange)./Q);
        quant_DCTblocks(irange,jrange) = quantblock;
    end
end
figure
imagesc(quant_DCTblocks)
colormap(gray)
set(gca,'Xtick',[])
set(gca,'Ytick',[])
axis('square')
title('quantized DCT Coeffs')

From the code it looks like the quantization step is doing two things: 1) dividing by "Q", and 2) rounding.  Undoing the divide by Q is easy.  The rounding is where the actual quantizing is happening though, and unfortunately that is the part that you cannot undo- at least not without more information.  Rounding causes information to be lost and you cannot undo that.


Image comparison using MATLAB I want to compare a bunch of depth images coming from Kinect and eliminate the similar ones. Those images present posture of a single user sitting on the chair (depicted in).

My idea is to create a function which will return the difference in range [0,1] between two postures. This means, I want to compare only a posture for given two images (a grey region). For example, if I pass images to my function, the result will be 0 (because postures are the same, but not at the same position on the picture).

Hence, I need to modify the following code (and find something to replace "correlation"):

a = imread('image_user_depth_000041.jpg');  %first img
b = imread('image_user_depth_000042.jpg');  %second img
a = RGB2IND(a,16);
b = RGB2IND(b,16);

c = corr2(a,b);           %finding the correlation
if c==1
 disp('The images are same') % display
else
 disp('the images are not same')
end;

Do you have any suggestions?
My regards and thanks!

You can borrow from techniques used by the image stitching community, where they have to deal with registering translated and rotated images. Your mistake was that you calculated the correlation before finding the optimal translation. The Wikipedia article should give you enough information to help yourself.

Mathematica Example

First, I split the picture above into frame1.png and frame2.png. Then I shifted the second one to align with the first, and trimmed both to retain their intersection. Finally, I compute the correlation coefficient by flattening the matrix into a vector.

img1 = ColorConvert[Import["frame1.png"], "Grayscale"]
img2 = ColorConvert[Import["frame2.png"], "Grayscale"]
im2t = ImageAlign[img1, img2, "Transformation" -> "Translation"]
imgtr = FindGeometricTransform[img1, img2, "Transformation" -> "Translation"]
{trx, try} = Part[TransformationMatrix@Last@imgtr, {1, 2}, 3]
{trim1, trim2} = ImageTake[#, {1, Floor@try}, {1, Floor@trx}] & /@ {img1, im2t}
{cimg1, cimg2} = Flatten@First@ImageData[#, "Interleaving" -> False] & /@ {trim1, trim2};
Correlation[cimg1, cimg2]

The result I got was 0.967834, indicating a close match. If necessary I would have considered scaling and rotation to better align the images but it was not necessary here.

Here is a MATLAB example.


Non-symmetrical DFT of real valued signal I have a homework assignment to take the $64$ DFT of $\cos((5\pi/32)n)$ and $\cos((5\pi/64)n)$. However, for the second case I obtained DFT spectrum such that I have that famous symmetry in real part and imaginary part of it because it is real valued signal. However, for the first signal I found something that has no symmetry in imaginary part of the DFT, but it's real too.

Why there is no symmetry although signal is real? Does it have anything to do with "windowing" the signal?

I wanted to calculate it analytically so that here is my code and the results:

....
t = 0:63;
x1 = cos(5*(pi/32)*t);
X1_analytic = zeros(1,64);
temp = 0;
for k = 0:63
    for n = 0:63
        temp = temp + x1(n+1)*exp(-1i*2*pi/64*k*n);
    end
    X1_analytic(k+1) = temp;
    temp = 0;
end
subplot(5,1,1);
stem(0:63,abs(X1_analytic),'--*r');
subplot(5,1,2);
stem(0:63,(X1_analytic+conj(X1_analytic))*(1/2));
title('Real part');
subplot(5,1,3);
stem(0:63,imag((X1_analytic)));
title('Imaginary part');
subplot(5,1,4);
stem(0:63,abs(X1_analytic));
title('Abs value ');
subplot(5,1,5);
stem(0:63,angle(X1_analytic));
title('Phase ');

Your choosen example frequency $\omega_0 = \frac{5 \pi} {32}$ touched a numerically sensitive case, in which the imaginary part of the theoretical DFT which should be zero is not observed so, on the numerical computation. Since the imaginary part is computed as non-zero, the phase which should be zero according to $\phi[k] = \tan^{-1}( \frac{\Im\{X_1[k]\}}{\Re\{X_1[k]\}} )$ whenever $\Im\{X_1[k]\}$ is zero, turns out to be some erratic function of $k$ (Note also that real part of $X_1[k]$ is also not computed to be exactly zero for those $k$ for which it should be, which also contributes to the erratic phase computation behaviour)

Consider the discrete time signal $x_1[n] = \cos(\frac{5\pi n}{32})$ and its 64-point DFT $X_1[k]$ (for ineteger $k \in [0:63]$) which is computed as:

$$ X_1[k] = \sum_{n=0}^{63} \cos(\frac{5\pi n}{32}) e^{-j \frac{2\pi}{64} nk} $$

expand the cosine into complex exponentials:
$$ X_1[k] = 0.5 \sum_{n=0}^{63} ({e^{j \frac{5\pi}{32} n} + e^{-j \frac{5\pi}{32} n}}) e^{-j \frac{2\pi}{64} nk} $$

adjust the cosine based complex exponential arguments:
$$ X_1[k] = 0.5 \sum_{n=0}^{63} ({e^{j \frac{10\pi}{64} n} + e^{-j \frac{10\pi}{64} n}}) e^{-j \frac{2\pi}{64} nk} $$

separate the sums:
$$ X_1[k] = 0.5 \sum_{n=0}^{63} {e^{j \frac{10\pi}{64} n} e^{-j \frac{2\pi}{64} nk} + 0.5 \sum_{n=0}^{63}e^{-j \frac{10\pi}{64} n}} e^{-j \frac{2\pi}{64} nk} $$

Now replacing numerical value of $N = 64$ with $N$ we shall see the summations:
$$ X_1[k] = 0.5 \sum_{n=0}^{N-1} {e^{j \frac{10\pi}{N} n} e^{-j \frac{2\pi}{N} nk} + 0.5 \sum_{n=0}^{N-1} e^{-j \frac{10\pi}{N} n}} e^{-j \frac{2\pi}{N} nk} $$

Merge the complex exponentials:
$$ X_1[k] = 0.5 \sum_{n=0}^{N-1} {e^{j \frac{2\pi n (5 -k)}{N}}+ 0.5 \sum_{n=0}^{N-1}e^{-j \frac{2\pi n (5+k)}{N}}} $$

Now the above sums add up to $N/2$ whenever $k=5$ or $k=-5$ ( this is also for $k = -5+ 64 = 59$) and they become zero for any other integer $k$ in the range $k \in [0:63]$

As you can now see, $X_1[k]$ is a real and even symmetric function of k. And its imaginary part is zero, however in your loop computation this
last condition is not realized due to numerical issues and therefore you observe erratic phase computation behaviour.


Closest two frequencies that can be resolved Forgive me if I misuse the terminology.

I wish to transmit two slowly changing analogue signals by amplitude modulating two higher frequency sinusoidal carrier waves. I will superimpose the two AM signals, send them down a wire, and then demodulate them at the other end using a digital signal processor.

To demodulate them, I will simply use two DFT bins.

I want to select two carrier frequencies as close together as possible. What is the closest I can get two carrier frequencies without them interfering?

There is no issue of frequency resolution in a typical modulation and demodulation scheme. (unless you want to implement some frequency/signal detectors)

Putting two message signals in adjacent channels requires that their spectrums do not overlap. In your case you should be able to specify the associated bandwidths of those two slowly changing analogue signals.

Assuming that associated (baseband) bandwidths of those signals are $F_1$ and $F_2$ Hz, then the carrier frequencies $f_{c1}$ and $f_{c2}$ have the relation:
$$ f_{c2} - f_{c1} > F_1 + F_2 $$

That's the minimum distance allowed for classical AM modulation, you could have used, instead, a single side band AM (SSB-AM) which would reduce the minimum distance to :
$$ f_{c2} - f_{c1} > \min(F_1,F_2)$$

But definitely that would require a more complex and difficult modulator and demodulator in terms of implementation.

Assuming that you have met this criterion (and assuming that you have selected the simple classical AM option), after down-converting, demodulating you should sample the signals satifying their respective nyquist rates.


Confusion in the derivation of convolution formula $$x(t)=\int_{-\infty}^{\infty}x(\tau)\delta(\tau-t)d\tau=\int_{-\infty}^{\infty}x(\tau)\delta(t-\tau)d\tau$$
Due to linearity,
$$y(t)=T\{x(t)\}=\int_{-\infty}^{\infty}x(\tau)T\{\delta(\tau-t)\}d\tau=\int_{-\infty}^{\infty}x(\tau)T\{\delta(t-\tau)\}d\tau$$
considering the impulse response of system to be h(t),due to time invariance of unit impulse signal ,
$$y(t)=\int_{-\infty}^{\infty}x(\tau)h(\tau-t)d\tau=\int_{-\infty}^{\infty}x(\tau)h(t-\tau)d\tau$$
here, i know i have done something wrong, because
$y(t)=\int_{-\infty}^{\infty}x(\tau)h(\tau-t)d\tau$
is the formula for cross-correlation between x(t) and h(t)
and $y(t)=\int_{-\infty}^{\infty}x(\tau)h(t-\tau)d\tau$ is the formula for convolution between x(t) and h(t)
but i cannot find, what exactly i have done wrong. could you help?

Your confusion stems from abusing the time invariance property.

Let's restate it for a simple 1D system:

$$ \text{if} ~~ T\{ x(t) \} = y(t) \implies T\{x(t-d) \} = y(t-d) $$ is correct but the following

$$ T\{ x(-t) \} = y(-t) ~~~\text{is wrong in general} $$

Now in your derivation you replace the input $\delta(t-\tau)$ with time reversed input $\delta(\tau-t)$ and claim the same on the output as $h(t-\tau) = h(\tau-t)$ which is not true as shown above.

The solution comes simply by recognizing the fact that $\delta(t-\tau) = \delta(-(t-\tau)) = \delta(\tau-t)$ and then proceed only with the $\delta(t-\tau)$ case as Marcus have stated.


Recognizing math functions within songs I'm new to DSP, and just discovered this StackExchange, so apologies if this isn't the right place to post this question.

Is there a resource that describes genres in a more mathematical terms? For example, if I've performed an FFT on the signal on this section of the song (2:09 if the link doesn't start there), is there any way of me being able to detect that this section has that rough sort of sound? Do sounds like this follow some mathematical function with which I can compare?
http://www.youtube.com/watch?v=SFu2DfPDGeU&feature=player_detailpage#t=130s (link starts playing sound straight away)

Is the only way to use supervised learning techniques, or is there a different approach (which  preferably doesn't require supervision)?

Thank you for any advice.

I think the distinction you're looking for is more like empirical vs. theoretical (as opposed to supervised vs. unsupervised), but I could be wrong about that.  In other words, the ideal thing would be to have a theoretical definition of various genres, rather than just a bunch of opaque data which can be used to classify a song [without any real understanding].
However, for general genre classification, you're probably stuck at least with training from examples, even if just to create the definitions of genres in the first place.  With respect to your example, consider how frequently people will argue [on YouTube] over whether a given track is really dubstep (e.g. any track that's more dubby and less wobbly, even though the genre started out without any real wobble). People define genres over time through examples, so it's reasonable to expect that algorithms which replicate that behavior would also require some examples.  The way people describe genres is almost like a feature vector anyway--they ask a list of questions about the song (e.g. Is it more breaky or wobbly? Does it have a lot of sub bass? How long is it? What's the tempo? Is there a vocal? etc.).
Of course, you may be able to choose a list of features that also provide an intuitive understanding of the genre.  A feature like "Dynamic Range" is something that a person can also detect by ear, but something like "Time Domain Zero Crossings" wouldn't be very intuitive--even if it works well for classification.  The following paper has quite a few features that might be interesting to you:

George Tzanetakis, Perry R. Cook: Musical genre classification of
audio signals. IEEE Transactions on Speech and Audio Processing 10(5):
293-302 (2002) link.

For measuring roughness, psychoacoustic roughness would be a good place to start, but it might not be sufficient to distinguish between dubstep leads and electro leads, for example.  For finer-grained distinctions, one thing to look into is timbre recognition.  The following thesis has a decent survey of techniques:

T. H. Park, “Towards automatic musical instrument timbre recognition,”
Ph.D. dissertation, Princeton University, NJ, 2004. link.

There's also a model related to perceptual roughness in Timbre, Tuning, Spectrum and Scale  which is used for constructing custom scales for arbitrary timbres.  The idea is that harmonics which are very close together produce beat frequencies which are perceived as dissonance.  Paraphrasing from Appendix F and E,

When $F$ is a spectrum with partials at frequencies $f_1,f_2,...,f_n$,
the intrinsic dissonance [assuming unit amplitudes] is
$$ D_F = 1/2 \space \sum_{i=1}^{n}{} \space \sum_{j=1}^{n}{\space d\left({|f_i - f_j| \over{\min(f_i,f_j)}} \right) } $$
where
$$d(x) = e^{-3.5 x} - e^{-5.75 x}$$
is a model of the Plomp-Levelt Curve.

It's used for measuring how pleasing a given chord is with respect to a timbre (by minimizing the dissonance).  I don't know if either roughness of the psychoacoustic variety, or intrinsic dissonance would be very fruitful for your purposes on their own, but they may be useful in combination with other metrics.
You'll probably have more luck classifying timbres mathematically than genres. For example, strings have even and odd harmonics, but a clarinet has only odd harmonics (cf. Sawtooth wave, Square wave). Dubstep wobble tends to be done with LFO-driven filters (low pass and/or formant filters), so something like Spectral Flux (see [Tzanetakis], above) might be a good starting point as a feature. However, I doubt anyone has studied mathematical classification of wobble yet ;)


Question regarding filter implementation - Audio EQ Cookbook Some of you may know the infamous Audio EQ Cookbook. You can take a look here.

We tried to implement the functions there to achieve a  implementation for a Lowpassfilter.
However, when we try to lower our variable cutoff-frequency, we have some issues. With decreasing cutoff-frequency the signal is increasing in volume. Reaching a distinct treshold at around 4000- 6000Hz the filter stops working, because it computes samplevalues that are ever increasing, reaching around 2000000 shortly after lowering the cutoff frequency.

We checked our code multiple times,and the functions from the Cookbook should be right. If someone wants to see the code i can hand it to you later.

For any hints or solution we are more than thankful.

EDIT: The important part of the code.

float Filter::processFilter(float sample)
{
   float unprocessed = sample;

    float processed = (unprocessed) + (computeB1() * x_n_1) +
        (computeB2() * x_n_2) - (computeA1() * y_n_1)
        -(computeA2() * y_n_2);

    x_n_2 = x_n_1;
    x_n_1 = unprocessed;

    y_n_2 = y_n_1;
    y_n_1 = processed;

    return processed;
}

With the values you provided ($Q=1$, $f_0=5355\,\text{Hz}$, and $F_s=48\,\text{kHz}$) I get the following filter coefficients from the Cookbook formulas:

b =

   0.11789
   0.23578
   0.11789

a =

   1.32248
  -1.52844
   0.67752

Note that the formulas give you denominator coefficients that are not normalized. Many implementations use $a_0=1$, which can be easily achieved by dividing numerator and denominator coefficients by $a_0=1.32248$.

The maximum pole radius of this filter is $0.72$, which means that we're not at all close to the unit circle, so stability shouldn't be an issue, no matter which implementation you use.

You can add more information if you like, but from what I know up to now, the most likely problem is a bug in your code.

Just test your filtering routine by using as input sequence x=1,0,0,...,, which should give the following output (the impulse response):

ans =

   8.9144e-02
   2.8131e-01
   3.6860e-01
   2.8188e-01
   1.3694e-01
   1.3860e-02
  -5.4140e-02
  -6.9672e-02
  -5.2786e-02
  -2.5313e-02
  -2.2122e-03
   1.0411e-02
   1.3166e-02
   9.8827e-03
   4.6767e-03
   3.4194e-04
  -2.0007e-03
  -2.4875e-03
  -1.8499e-03
  -8.6361e-04
  -5.0389e-05


What is the minimum sampling rate for QPSK or what is the bandwidth of QPSK signal? Consider  we have created a QPSK signal with a carrier frequency 100kHz. Before applying matched filter or demodulation we need to take samples from the received signal. So what is the minimum sampling rate we must have? Consider that somehow we have perfect synchronization and there is no drift of the carrier.

Any help would be appreciated.

I have asked this question (for more attention) at stackexchange .

As best as I can tell from the terminology, QPSK a.k.a. 4-PSK or 4-QAM can convey twice the number of bits of information as BPSK in the same bandwidth.

Given a one-sided bandwidth of $B$ and assuming sufficient S/N ratio so that no bits are corrupted, BPSK can transfer bits at a rate of $2B$.  That means, given the same one-sided bandwidth, QPSK can transfer bits at a rate of $4B$.

But the sample rate must be at least $2B$.  If the receiver is perfectly synchronized to the transmitter, there would be one sample per bit with no bit error, if the transmitter sends a bipolar $\operatorname{sinc}(\cdot)$ for each bit.

But for QPSK, each sample is a complex value with a real part and imaginary part, each part containing one bit sample if this is 4-PSK.

So if you use quadrature demodulation to bump the IF or RF down to baseband $i(t)$ and $q(t)$, and sample those signals, you will need one complex-valued sample per pair of bits.


Why is a circular mask appropriate for Fourier filtering rectangular images? Suppose I apply 2D DFT to an image with dimensions $H{\times}W$ where $H \neq W$, then shift the DC component to the center. Why does a circular mask capture the lowest frequency components, i.e. why is it not an ellipse given that the image is rectangular? My concern is that for rectangular images, the K lowest frequencies might be arranged in a non-circular pattern.

For simplicity, let's not do any shifting and only consider non-negative frequencies.
Let's assume that the horizontal and vertical image dimensions are even integers $W$ and $H$. Looking at the output of a $H\times W$ 2-d DFT of the image, the $u$th column with $u\le W/2$ represents a horizontal frequency of $u/W$ times the horizontal sampling frequency and the $v$th row with $v\le H/2$ represents a vertical frequency of $v/H$ times the vertical sampling frequency. For a square image grid the horizontal and vertical sampling frequencies are equal and in the following denoted by a single variable $f_s$. The frequency-magnitude of a 2-d frequency at bin $u, v$ will then be $\sqrt{(u/W)^2 + (v/H)^2}f_s$.
For a cut-off frequency $f_c$ your mask would select frequencies:
$$\sqrt{\left(\frac{u}{W}\right)^2 + \left(\frac{v}{H}\right)^2}f_s < f_c\tag{1}$$
$$\Rightarrow\frac{f_s^2}{W^2f_c^2}u^2 + \frac{f_s^2}{H^2f_c^2}v^2 < 1.\tag{2}$$
That indeed defines an ellipse in coordinates $u, v$.
However, if you consider the actual frequencies $\frac{u}{W}f_s$, $\frac{v}{H}f_s$ as coordinates, then what you have is a circle:
$$\text{Eq. 1}$$
$$\Rightarrow\frac{1}{f_c^2}\left(\frac{u}{W}f_s\right)^2 + \frac{1}{f_c^2}\left(\frac{v}{H}f_s\right)^2 < 1.\tag{3}$$
To summarize, it depends on how you express your frequencies.


Processing OBPSK as OQPSK I'm trying to create and process an OBPSK signal, which, as I understand it, transmits even bits on one axis and the odd bits on the opposite axis, resulting in a constellation with 4 unique points. OBPSK seems to be somewhat rare, and a lot of software doesn't have functions to handle it. I created it (if I understand it correctly) by modulating the even and odd bits as BPSK separately, and rotating one of them by 90 degrees. However, with 4 constellation points, it seems similar to QPSK, or more specifically OQPSK, which is common with software having built in functions for modulation, etc. Is there some way to process OPSK, perhaps with some sort of transformation and/or change in data rate, as OQPSK?

OBPSK in the way that you describe it is also called $\pi/2$-BPSK meaning that the signal constellation is different from bit to bit, rotating by $\pi/2$ after it is used. There are two variations of this format.

Use $\pm \cos(\omega_0 t)$ to transmit one bit and
$\mp \sin(\omega_0 t)$ to transmit the next bit. Repeat this over and over again.  Note the $\mp$ before the $\sin$ carrier.
This method rotates the constellation by $\pi/2$ counterclockwise in going from $\pm \cos(\omega_0 t)$ to $\mp \sin(\omega_0 t)$, and clockwise when going from $\mp \sin(\omega_0 t)$ to $\pm \cos(\omega_0 t)$.
Use

$\pm \cos(\omega_0 t)$ to transmit the first bit
$\mp \sin(\omega_0 t)$ to transmit the second bit
$\mp \cos(\omega_0 t)$ to transmit the third bit
$\pm \sin(\omega_0 t)$ to transmit the fourth bit,
and repeat this over and over. Note the $\mp$ and $\pm$ signs. This
corresponds to the constellation rotating by $\pi/2$ counterclockwise after each and every bit transmission.

Neither form of OBPSK signal can be demodulated very well
using a OQPSK demodulator, though it would be possible to do
so. In OQPSK,
even-numbered bits are transmitted using $\pm\cos(\omega_0t)$ and
odd-numbered bits using $\mp\sin(\omega_0 t)$ but the important
difference is that the $2n$-th bit is transmitted during the
time interval
$[2nT, (2n+2)T)$ while the $(2n+1)$-th bit is transmitted during
the time interval $[(2n+1)T, (2n+3)T)$, i.e. in overlapping but offset time intervals.
At the demodulator, the matched filters/correlators in the I branch
and Q branch thus integrate/correlate over offset time intervals,
and the demodulated data bits come out of the OQPSK demodulator alternately on the I and Q lines, one every T seconds.
On the other hand, in
OBPSK, the $2n$-th bit is transmitted during $[2nT, (2n+1)T)$ and
the $(2n+1)$-th bit during $[(2n+1)T, (2n+2)T)$. If demodulated with
a plain vanilla OQPSK demodulator, during $[(2n+1)T, (2n+2)T)$ the matched filter/correlator input in the I branch would be just noise
(and similarly for the Q branch).
This will decrease the effective SNR. So, Yes, a standard OQPSK
receiver could be used to demodulate an OBPSK signal, but
there are better ways of demodulating an OBPSK signal than just
processing it through an OQPSK receiver.