Real discrete Fourier transform I am trying to understand the real DFT and the DFT and why the distinction exists.

From what I know so far the DFT uses $e^{i2\pi kn/N}$ for basis vectors and gives the representation $$x[n]=\sum_{k=0}^{N-1}X[k]e^{i2\pi kn/N}$$ The sum is written from $k=0$ to $N-1$ for historical reasons I think instead of writing it in a way analogous to the Fourier series with the sum going from $k=-N/2$ to $N/2-1$: $$x[n]=\sum_{k=-N/2}^{N/2-1}X[k]e^{i2\pi kn/N}$$
This relying on a peculiar anomoly of the DFT where high frequencies are the same as negative frequencies: $e^{i2\pi kn/N}=e^{i2\pi (k-N)n/N}$.

Continuing the analogy with Fourier Series the real DFT gives the representation
$$x[n]=\sum_{k=0}^{N/2}\left(X_R[k]\cos\left(\frac{2\pi kn}{N}\right)-X_I[k]\sin\left(\frac{2\pi kn}{N}\right)\right)$$
This can be viewed as pairing $e^{i2\pi kn/N}$ with $e^{-i2\pi kn/N}$ in the DFT representation where the sum ranges from $k=-N/2$ to $N/2-1$. This is very much like the pairing $c_n e^{in\theta}+c_{-n}e^{-in\theta}=a_n \cos n\theta + b_n \sin n\theta$ which connects the two representations of a Fourier Series:$$\sum_{-\infty}^\infty c_n e^{in\theta}= \frac{a_0}{2} + \sum_1^{\infty}(a_n \cos n\theta + b_n \sin n\theta)$$

My question then is why is the DFT so much more prevalent than the real DFT? One would expect that since the real DFT is using real valued sines and cosines as the basis and is thus representing the geometric picture better that people would like it more. I can see why the DFT and the continuous Fourier Transform would be preferred in a theoretical sense as the algebra of exponentials is simpler. But ignoring the simpler algebra, from a practical computational applied viewpoint why would the DFT be more useful? Why would representing your signal with complex exponentials be more useful in various physics, speech, image, etc. applications than decomposing your signal into sines and cosines. Also if there is anything subtle I'm missing in my above exposition I would like to know: I'm puzzled that the DFT is seen as more connected to the continuous Fourier Transform than to the Fourier Series.

The advantage of the complex DFT or complex Fourier transform or
complex Fourier series is that linear
systems have the nice property that the response to $A\exp(j\omega t)$ is
$H(\omega)A\exp(j\omega t)$.  (Here $A$ can be a complex constant).
So the output is just a scalar multiple of the input.
More importantly, if we have a representation of the input
as a weighted sum of complex exponentials, the output is
just another weighted sum of the same exponentials.  Different
weights, but same set of exponentials. Furthermore, each new weight
is obtained by multiplying the old weight by an appropriate number.

Of course, no physical system has complex-values
signals going in and coming out; at least, not as of today though
one can always hope for better things in the future.  In the mean time,
we take real parts of the complex signals, or get the response to $\cos(\omega t)$ or
$\sin(\omega t)$ via linearity and superposition and liberal use of
$$\begin{align*}
\cos(\omega t) &= \frac{\exp(j\omega t) + \exp(-j\omega t)}{2}\\
\sin(\omega t) &= \frac{\exp(j\omega t) - \exp(-j\omega t)}{2j}
\end{align*}$$

In contrast, the response to $\cos(\omega t)$ is of the form
$B(\omega)\cos(\omega t) + C(\omega)\sin(\omega t)$.  So,
while linearity and superposition etc all work, the output
might well need the use of different basis functions than the
input does.  Very closely related, of course, but still possibly
different and maybe more  basis functions might be needed. For
example, input $\cos(\omega t)$ is represented by one basis function,
output $B(\omega)\cos(\omega t) + C(\omega)\sin(\omega t)$ by
two basis functions.  It can be argued that complex functions
require twice as much work as real functions and so any savings
are purely imaginary (pun intended), but complex representations allow
uniform treatment while sin/cos representations do not.  Quick!
Given the response to $\cos(\omega t)$ is
$B(\omega)\cos(\omega t) + C(\omega)\sin(\omega t)$,
what is the response to $\sin(\omega t)$?  You have to work at it
a bit, you may need to invoke formulas such as
$$\cos (\alpha + \beta) = \cos(\alpha)\cos(\beta) - \sin(\alpha)\sin(\beta)$$
and so on.  With complex exponentials, life is a lot easier.

But, as in real life, your mileage may vary, and if you feel
that sin/cos representations are the way to go and complex
exponentials should be eschewed, you are free to follow your heart.
If you have difficulty communicating your ideas to colleagues,
bosses, clients or consultants, that will be their loss, not yours.


Figuring out how an RF pulse affects the MRI spins So, I understand how the slice selection gradients work in MRI. So, the frequency offset introduced by the slice selection gradient at a location $z$ relative to the MRI isocenter is given by:

$$
\Delta f = \gamma z G_z
$$

Putting some numbers to it, for a gradient of 5 mT/m and a slice thickness of 3 mm, we introduce an offset of $639.7$ Hz per slice.

Now, say I have an RF pulse which has a bandwidth of 2000 Hz and a center frequency of 10 kHz. Given the location of the spin along the slice direction, I want to know whether this spin will be excited by the RF pulse or not.

Is it possible to figure that out from this information?

I have a feeling it should be possible from the relation above but I am not sure how the bandwidth parameter will come into play here.

Is it a simple matter of checking that if the resonant frequency experienced by the spin falls in the RF center frequency $\pm$ bandwidth / 2.0?

Yes, it basically is that simple. If the resonance frequency of the spin falls in the frequency range that the pulse excites, the spin will be manipulated by the RF pulse. Reality usually kicks in when you have imperfect slice profiles (i.e. to achieve a perfect rectangular slice profile, you would need an infinitely long RF sinc-pulse).

Is it possible that you accidentally mixed up angular frequency and normal frequency? I am asking, because your frequency spread and RF bandwidth would excite $\frac{2000\,\textrm{Hz}}{639.7\,\frac{\textrm{Hz}}{\textrm{slice}}} \approx 3.12\,\textrm{slices}$, if I am not mistaken, and... when it comes to frequencies and a ratio is so close to $\pi$, and there is a "/2.0" in the original question, there could be something going on with wrong frequency units.


Fourier transform of a signal and its autocovariance function While I do know the difference between the two, in theory, I am not very sure about why we look at the Fourier transform of the auto-covariance function. What extra information does it give us over and above what we can already get by taking the Fourier transform of the original time-series? Also, when is one advantageous over the other?

The interest in looking at the Fourier transform of the autocorrelation (which is usually equivalent to autocoveriance) comes from the Wiener-Khinchin theorem that states
Consider a wide-sense stationary signal $x(t)$. Then the power spectral density $S_x(f)$ is given by
$$S_x(f) = \int_{-\infty}^{\infty}R(\tau)e^{-j2{\pi}f{\tau}}d{\tau}$$
Where $R(\tau)$ is the autocorrelation function of $x(t)$.
There are times where the Fourier transform of a signal does not exist. However, this result still allows us to analyze the signal in the frequency domain and retrieve information by using the autocorrelation and not the signal itself.
This is also helpful when analyzing LTI systems. Given a system $h(t)$ with input $x(t)$ and output $y(t)$, we have a relationship between the power spectral densities given by
$$S_y(f) = S_x|H(f)|^2$$
This aids in estimating properties of the signals as well as the system itself, depending on what information is already known.
In the discrete domain, the DFT can always be calculated for a finite-length signal. The theorem still holds however, and despite being able to "take" the Fourier transform as in the DFT, errors may begin to appear when estimating parameters.
References by OP's request:
https://en.wikipedia.org/wiki/Wiener%E2%80%93Khinchin_theorem
https://mathworld.wolfram.com/Wiener-KhinchinTheorem.html
https://ocw.mit.edu/courses/electrical-engineering-and-computer-science/6-011-introduction-to-communication-control-and-signal-processing-spring-2010/readings/MIT6_011S10_chap10.pdf


State-of-art of telecommunications? There seems to be a number of different ways for telecommunication. What is the most recent (and maybe standard) technique used today? It seems the simple AM/FM modulations are definitely somewhere in the past, then there are these spread spectrum systems... All hints welcome :)

You should realize that when you say "telecommunications" you are talking about a huge field involving lots of different technologies.  Based on the mention of AM/FM modulations I infer that you are referring to digital (though AM/FM are actually analog) wireless transmission techniques (though I may have inferred incorrectly since wired modems can use some pretty sophisticated modulations).  Even restricting ourselves to wireless telecommunications, there are lots of different technologies that we could talk about- modulation types are just one.

Anyway, back to what I think you were asking.

The progression of modulation types went something like the following:

AM (very simple analog)
FM (more sophisticated analog- but uses more bandwidth than AM)
OOK (very simple digital)
FSK (digital version of FM)
PSK
QAM
TCM
DSSS
OFDM

If you google any of the above terms and add "modulation", you should get some useful references.

OFDM is, as far as I am aware, the latest-and-greatest modulation type.  It has the useful feature of being extremely efficient at packing data into a given amount of bandwidth.  Since bandwidth is at a premium for the telecoms, that is very important.  Many of the other modulation types are still used though, and will continue to be used because they have other useful features.


A 10-seconds digital reverb with very few memory (64KB) I'm currently thinking about the design of a low-cost digital reverb guitar pedal, using a Teensy or Arduino.

My goal is to achieve big 10-seconds long reverb.

If I implement reverb with a mix of many delay lines, this will require to store in RAM at least 10*44100*2bytes ~ 0.8 MB of audio data (44.1Khz, 16bits, mono), so it seems to be impossible to achieve with a 64KB RAM board (am I wrong?).

My question: Is there a more clever way (using FFT ? overlap-add stuff, or anything else?) to implement 10 seconds reverb with very little RAM (64KB)?

Take a look at this https://ccrma.stanford.edu/~jos/pasp/Freeverb.html or https://github.com/highfidelity/gverb. Gverb is clearly the better sounding one but harder to implement.

Not sure if you can squeeze it into 64kB but it's worth a shot. The length of the delay lines is scalable to some extent. It roughly corresponds to room size with small delays being small room. 10s in a small room sounds more like bathroom and less like a concert hall or a cathedral. For really high reverb times I would recommend a floating point implementation to avoid instability and noise in the feedback loops.


Find original Hz of converted audio-file? I have received an audio file which is 44.1kHz wav. My ears suspect this audio clip was originally lower than 44.1, and encoded to 44.1 later, much like re-encoding a 720p video to a 1080p video. Is there a way to find out if this is true? Like, analyse the audio clip to see.. how many Hz it uses..? I have literally no idea how this works.

The short answer is no.

The longer answer is that if there is some known signal that you know the frequency you could but this is knowledge you bring to the problem.

You can look around 60Hz or 50Hz depending what the utility power frequency uses and see it there is some suggestion of a steady line in the data in that frequency range.

Musicians, after thousands of hours of playing an instrument will know what it should sound like.  While somewhat subjective, musical aesthetics are more than adequate for aesthetic purposes.

Summary,
If all you have is a file, not a chance.

If you have some clues and hunches and familiarity with the data there is a reasonable possibility.


How to calculate the Fourier transform of absolute-valued sinusoid? The signal $x(t) = |1+a\sin(\omega t)|,(a>0)$ is a continuous waveform. In order to extract the frequency parameter $\omega$, I conduct the FFT of it and obtain its spectrum showed as follows. The frequency $\omega$ is set as 1 Hz. Spikes occur at the fundamental frequency and its harmonics. I'm curious that why this phenomenon exists. How to calculte the Fourier transform of $x(t)$?

Spectrum of $X(f)$

The waveform of $x(t)$

For $0 < a \le 1$, the Fourier Transform is straightforward:

$$\begin{align*}\mathscr{F}\left\{|1+a\sin(2\pi s_0t)|\right\} &= \mathscr{F}\left\{1+a\sin(2\pi s_0t)\right\} \\
&= \int_{-\infty}^{\infty}\left[1+a\sin(2\pi s_0t)\right]e^{-2\pi i st} \mathrm{d}t\\
& = \delta(s) +\dfrac{a}{2i}\left[\delta(s-s_0)-\delta(s+s_0)\right] \\
\end{align*}$$

For $ a > 1$, we have to rewrite the function to get rid of the absolute value bars:

$$|1+a\sin(2\pi s_0t)| = |s_0|\mathrm{III}(s_0t) *\left(\left[1+a\sin(2\pi s_0t)\right]\cdot\left[\mathrm{\Pi}\left(\dfrac{\pi s_0}{\frac{\pi}{2}-\mathrm{Arcsin}\left[-\frac1{a}\right]}t-\dfrac{\frac{\pi}{2}}{2\pi s_0}\right)-\mathrm{\Pi}\left(\dfrac{\pi s_0}{\frac{\pi}{2}+\mathrm{Arcsin}\left[-\frac1{a}\right]}t+\dfrac{\frac{\pi}{2}}{2\pi s_0}\right)\right]\right)$$

Where $*$ denotes convolution, $\mathrm{Arcsin}()$ denotes the principal arcsine,

$$\mathrm{III}(at) = \dfrac1{|a|}\sum_{n=-\infty}^{\infty}\delta\left(t-\dfrac{n}{a}\right)$$
and
$$\begin{equation}\mathrm{\Pi}(at) = \begin{cases}
1 & |at| <\frac1{2}\\
0 & |at| >\frac1{2}\end{cases}\end{equation}$$

Note that the $\mathrm{\Pi}()$ functions are two rectangles snipping out the two different segments of one cycle, and the convolution with the $\mathrm{III}()$ function creates all the periodic replicas of those two segments of that one cycle.

Now using tables of Fourier Transforms, and various, basic Fourier Transform theorems we get, for $a >1$,

$$\mathscr{F}\left\{|1+a\sin(2\pi s_0t)|\right\} =\mathrm{III}\left(\dfrac{s}{s_0}\right)\cdot\left(\left[\delta(s)+\dfrac{a}{2i}\left(\delta(s-s_0)-\delta(s+s_0)\right)\right]*\left[\dfrac{\frac{\pi}{2}-\mathrm{Arcsin}\left(-\frac1{a}\right)}{\pi|s_0|}e^{2\pi i\left(\dfrac{-\frac{\pi}{2}}{2\pi s_0}\right)\left(\dfrac{\frac{\pi}{2}-\mathrm{Arcsin}\left(-\frac1{a}\right)}{\pi s_0}\right)s}\mathrm{sinc}\left(\dfrac{\frac{\pi}{2}-\mathrm{Arcsin}\left(-\frac1{a}\right)}{\pi s_0}s\right)-\dfrac{\frac{\pi}{2}+\mathrm{Arcsin}\left(-\frac1{a}\right)}{\pi|s_0|}e^{2\pi i\left(\dfrac{\frac{\pi}{2}}{2\pi s_0}\right)\left(\dfrac{\frac{\pi}{2}+\mathrm{Arcsin}\left(-\frac1{a}\right)}{\pi s_0}\right)s}\mathrm{sinc}\left(\dfrac{\frac{\pi}{2}+\mathrm{Arcsin}\left(-\frac1{a}\right)}{\pi s_0}s\right)\right]\right)$$

So you have the difference of two modulated $\mathrm{sinc}()$ functions, with 1 real replica at DC and 2 imaginary replicas at $\pm s_0$ (due to convolution with the $\delta()$ functions), sampled at integer multiples of $s_0$ (by the $\mathrm{III}()$ function).


Can some one give me an example of cross correlation's high value where template matches with image? I find it difficult to understand how result of cross correlation of a patch of image and image is high at matching point? what if patch contains small values and other non matching point contains very high values.

Please read this document here. If I understood your question clearly, then you need to take out the mean from the image. After that probably you would find cross-correlation for the matching patch.


What does it mean for an image to be "Markovian"? I've been following the recent lectures on Coursera for image processing. In the lecture on the use of the discrete cosine transform  (DCT) for image compression, he says that for images that are "Markovian", the DCT has the same properties as the Karhunen-Loeve Transform. He uses this argument to support the use of the DCT in general image compression, since small subsets of natural images usually have this property. What does "Markovian" mean in this context? Is he referring to Markov Random Fields where the graph is constructed by pixels as nodes and edges between neighboring pixels? What are some examples of images that violate this property?

It's a probabilistic model that is assumed to generate the image. To make theoretical computations on the performance of certain image processing tasks such as DCT transforming the image for data compression, the image is modelled with a simple mathematical equation.

The Markov image generation model, is used to produce images where there will be high enough correlation between neighbouring samples. This high correlation is therefore efficiently utilized by a transform such as DCT to achieve compression ratios that is close to the maximum available with a nonpractical KLT (Karhunen Loeve Transform)

As far as I know, it is not much related with markovian fields.
The simplest 1-D example of such a markov source model, assuming linear dependence between samples, is like: $$x[n] = \rho x[n-1] + {\omega}[n]$$

This simple dependency clearly shows the fact that current sample x[n] is produced from past sample x[n-1] plus a white noise sample w[n], where $\rho$ is the correlating factor and as it is close to 1, samples are highly correlated, when it is close to zero, samples become more noisy and uncorrelated.

For a real image, such a model can only predict samples in a very short duration, because of the non-stationarity of the image. This is one of the reasons why DCT based transforms are applied for short lengths of blocks such as 8x8 pixels.


Deterministic / Non-deterministic Stochastic Process Problem 6.1-6 of Probability, Random Variables, and Random Signal Principles, 4th Edition by Peebles asks

  If a  process is defined by $X(t) = A$, where $A$ is a continuous random variable uniformly distributed on  $(0,1)$,   is
  the process deterministic or non-deterministic?

The textbook says that  "a process is called deterministic if its value as a function of time can be pre-determined" and the answer to the
question above is the process is deterministic.

My questions are

If $A$ is a random variable, then how can we determine time functions (or ensemble members) before specific $A$ is given?
What are examples of non-deterministic processes?

Regardless of what your book might say, a random process is a collection of random variables that is described in mathematical notation as
$\{X(t) \colon t \in \mathbb T\}$ or $\{X_t\colon t \in \mathbb T\}$ where
$\mathbb T$ is called the index set (typically, $\mathbb T$ is
$(-\infty,\infty)$ or $\{\cdots, -2, -1, 0, 1, 2, \cdots\}$ or subsets
thereof) and we have one random variable for each $t \in \mathbb T$.
Each random variable is, of course, a mapping from the sample space
$\Omega\to \mathbb R$, and when the experiment is performed and
the outcome $\omega \in \Omega$ is known, the random variable
$X_t$ maps that $\omega$ onto some real number that is, of
course, usually denoted by $X_t(\omega)$.

Now that we have this (typically infinite) collection of random variables,
a sample path or realization of the random process is the "waveform"
that you see when the experiment is performed (with outcome
$\omega_{11809}$, say)
and each of the random
variables $X_t$ maps this outcome onto $X_t(\omega_{11809})$.
That's the "waveform": its value at $t = 0$ is $X_0(\omega_{11809})$,
and when $t = 12$, its value is $X_{12}(\omega_{11809})$ and so on.

Now we come to the crux of the matter. A deterministic random process
is one for which if we know the value(s) of one or more of the
$X_t(\omega)$, say $X_{0}(\omega)$ and $X_2(\omega)$,
then we know the value of all the $X_t(\omega)$: knowing the
value taken on by a few of the random variables in the set tells
us the values that all the random variables took on. That is,
there is "no more randomness left" in the process.

With this in mind, consider a random process for which all the
random variables $X_t$ are in fact the same random variable
$A$. (Whether $A$ is uniformly distributed on $(0,1)$ or not is
totally irrelevant to the issue being considered). Then, if we
know that $X_0(\omega) = A(\omega) = 0.13$, then
we know that $X_t(\omega) = 0.13$ for all $t \in \mathbb T$.
That is why your book is claiming that the random process
is deterministic. Note that if $X_t$ were to equal
$A \cos(2\pi f_0 t)$ where$f_0$ is some known constant instead of just $A$,
the process would still be deterministic since knowledge of the
value of $X_0$ would give away knowledge of everything.
If instead $X_t = A \sin(2\pi f_0 t)$, then knowing the value
of $X_0$ would not be useful but knowing the value of
$X_{(4f_0)^{-1}}$ (or any other $X_t$ where $\sin(2\pi f_0 t) \neq 0$)
would again be a telltale.

As a harder question that you can use to test whether you understand
the above explanation, consider whether a random process in which
by $$X_t = A \cos(2\pi f_0 t + \Theta), -\infty < t < \infty$$ where $A$ and $\Theta$
are independent random variables is a deterministic process or not.
Does knowing the values of just a few of the random variables $X_t$
tell you the value of all of them?

So, that is the basic idea behind the notion of deterministic
random processes: observation of some part of a sample path
allows us to determine the entire sample path with perfect
accuracy. How about nondeterministic random processes? Consider
$$X_n = A_n, \quad n = \cdots, -2, -1, 0, 1, 2, \cdots$$
where the $A_n$ are independent random variables with some
specific distribution, e.g. each $A_n \sim \mathcal N(0,\sigma^2)$.
In this case, the process is often referred to as a
(discrete-time) white noise process which might give you a
hint as to whether it should qualify as a deterministic or
a nondeterministic process.

Exercise: what is the significance of the number 11809 in the second
paragraph above?


Finding the $\mathcal Z$-transform of $((-\frac{1}{3})^n + \frac{1}{2})^n \mu[n-2]$? What is the $\mathcal Z$-transform of $\left(\left(-\frac{1}{3}\right)^n + \frac{1}{2}\right)^n \mu[n-2]$?

Doing raw computations with large $n$ gives a sum of $0.6030$ which doesn't seem right, and I'm not sure how to manipulate the expression to get some clean result.

Since $\left| -\frac{1}{3}\right|< \frac{1}{2}$ and $n$ starts at $2$, one can think about an approximation based on the first terms of the polynomial development of $a_n = \left(\frac{1}{2}\right)^n\left(1+2\left(-\frac{1}{3}\right)^n\right)^n $. Using a binomial (polynomial) development up to degree $2$ on the second term only, one gets:

$$a_n \approx b_n= \left(\frac{1}{2}\right)^n + 2n\left(-\frac{1}{6}\right)^n + 2n(n-1)\left(\frac{1}{18}\right)^n \,.$$

Details: since coefficients are zero for $n<2$, for $n=2$,

$$a_2=\left(1+2\left(-\frac{1}{3}\right)^2\right)^2=  b_2 = 1+ {2 \choose 1}\left(2\left(-\frac{1}{3}\right)\right)+ {2 \choose 2}\left(2\left(-\frac{1}{3}\right)\right)^2$$

exactly. Then for $n>2$, the term $c_n=2\left(-\frac{1}{3}\right)^n$ tends to be small. You can approximate $(1+c_n)^n$ by  $1+ {n \choose 1}c_n+ {n \choose 2}c_n^2$.

The following graphs provide an evaluation of the approximation, up to $n=8$. Top graph superimposes $a_n$,  @MaximilianMatthé's approximation (MM), and what we get using the first (D0), the first two (D1), or the three terms (D2) of the above series. $D0$ (only $\left(\frac{1}{2}\right)^n$) is clearly insufficient. The others are quite close. The bottom graph shows the absolute approximation errors. As you can see, $D1$ and $D2$ converge very fast when $n$ grows, and are already quite good for the first terms.

The nice thing with these approximations is that you not only reuse the classical properties of  the $\mathcal{Z}$-transform: linearity, time-shift and scaling, but also differentiation (up to order $2$), because of the factors $n$ and $n(n-1)$:

$$ n x_n \to -z\frac{dX(z)}{dz}\,.$$

A Matlab code for reproduction:

n = (2:8)';
s = ((-1/3).^n+1/2).^n;
s0 = ((-1/3).^n+(1/2).^n);
sa = ((1/2).^n);
sb = ((1/2).^n).*(1+2.*n.*(-1/3).^n);
sc = (1/2).^n+2.*n.*(-1/6).^n + 2.*n.*(n-1).*((1/18).^n);

subplot(2,1,1)
plot(n,[s,s0,sa,sb,sc],'.-');axis([0 max(n) -Inf Inf]); grid on
legend('Orig.','MM','D0','D1','D2')
h=xlabel('$n$ index');set(h,'INterpreter','latex');title('Series')
subplot(2,1,2)
semilogy(n,abs([nan(size(s)),s0-s,sa-s,sb-s,sc-s]),'.-');axis([0 max(n) -Inf Inf]); grid on
legend('Orig.','MM','D0','D1','D2')
h=xlabel('$n$ index');set(h,'INterpreter','latex');title('Absolute approximation errors')


Time domain maximum from frequency domain data? Is it possible to calculate the maximum value of a time-domain signal from frequency-domain representation without performing an inverse transform?

Suppose that Alice has a vector $\mathrm x \in \mathbb R^n$. She computes the DFT of $\mathrm x$

$$\mathrm y := \mathrm F \mathrm x \in \mathbb C^n$$

where $\mathrm F \in \mathbb C^{n \times n}$ is a Fourier matrix. Alice then tells Bob what $\mathrm y$ is. Since the inverse of the Fourier matrix is $\mathrm F^{-1} = \frac 1n \, \mathrm F^*$, Bob can recover $\mathrm x$ via

$$\mathrm x = \frac 1n \, \mathrm F^* \mathrm y$$

and then compute $\| \mathrm x \|_{\infty}$ to find the maximum absolute value of the entries of $\mathrm x$. What if computing matrix inverses and Hermitian transposes is not allowed? Bob can then write $\mathrm F$ and $\mathrm y$ as follows

$$\mathrm F = \mathrm F_{\text{re}} + i \,\mathrm F_{\text{im}} \qquad \qquad \qquad \mathrm y = \mathrm y_{\text{re}} + i \,\mathrm y_{\text{im}}$$

and, since $\mathrm x \in \mathbb R^n$, the equation $\mathrm F \mathrm x = \mathrm y$ yields two equations over the reals, namely, $\mathrm F_{\text{re}} \, \mathrm x = \mathrm y_{\text{re}}$ and $\mathrm F_{\text{im}} \, \mathrm x = \mathrm y_{\text{im}}$. Bob can then solve the following linear program in $t \in \mathbb R$ and $\mathrm x \in \mathbb R^n$

$$\begin{array}{ll} \text{minimize} & t\\ \text{subject to} & - t 1_n\leq \mathrm x \leq t 1_n\\ & \begin{bmatrix} \mathrm F_{\text{re}}\\ \mathrm F_{\text{im}}\end{bmatrix} \mathrm x = \begin{bmatrix} \mathrm y_{\text{re}}\\ \mathrm y_{\text{im}}\end{bmatrix}\end{array}$$

which can be rewritten as follows

$$\begin{array}{ll} \text{minimize} & \begin{bmatrix} 1\\ 0_n\end{bmatrix}^{\top} \begin{bmatrix} t\\ \mathrm x \end{bmatrix}\\ \text{subject to} & \begin{bmatrix} -1_n & \mathrm I_n\\  -1_n & -\mathrm I_n\end{bmatrix} \begin{bmatrix} t\\ \mathrm x \end{bmatrix} \leq \begin{bmatrix} 0_n\\ 0_n\end{bmatrix}\\ & \begin{bmatrix} 0_n & \mathrm F_{\text{re}}\\ 0_n & \mathrm F_{\text{im}}\end{bmatrix} \begin{bmatrix} t\\ \mathrm x \end{bmatrix} = \begin{bmatrix} \mathrm y_{\text{re}}\\ \mathrm y_{\text{im}}\end{bmatrix}\end{array}$$

and not only recover $\mathrm x$ but also obtain $t = \| \mathrm x \|_{\infty}$. However, is solving a linear program cheaper than computing a Hermitian transpose?

MATLAB code

The following MATLAB script

n = 8;

% build n x n Fourier matrix
F = dftmtx(n);

% -----
% Alice
% -----

% build vector x
x = randn(n,1);

% compute DFT of x
y = F * x;

% ---
% Bob
% ---

% solve linear program
c = eye(n+1,1);
A_in = [-ones(n,1), eye(n);
        -ones(n,1),-eye(n)];
b_in = zeros(2*n,1);
A_eq = [zeros(n,1), real(F);
        zeros(n,1), imag(F)];
b_eq = [real(y); imag(y)];
solution = linprog(c, A_in, b_in, A_eq, b_eq);

% extract t and x
t = solution(1);
x_rec = solution(2:n+1);

% check results
disp('t = '); disp(t);
disp('Infinity norm of x = '); disp(norm(x,inf));
disp('Reconstruction error = '); disp(x_rec - x);

produces the output

Optimization terminated.
t =
2.2023

Infinity norm of x =
2.2023

Reconstruction error =
1.0e-013 *

0.0910
0.0711
0.0167
-0.1077
0.1049
0.0322
0.1130
0.2776

The original vector is

>> x

x =

-1.1878
-2.2023
0.9863
-0.5186
0.3274
0.2341
0.0215
-1.0039


Constant Modulus Algorithm and the Gradient Operation CMA is a blind channel equalization algorithm with the details presented above. I am facing 3 difficulties and shall appreciate help

Q1: Does $H$ and the bar over $\bar{y_k}$ represent the Transpose symbol?

Q2: When calculating the gradient, $\Delta J(w) = 2E{({|y_k|}^2 - 1). \Delta(w^Hx_kx_k^Hw)}$, how does the gradient $\Delta(w^Hx_kx_k^Hw)$ becomes $x_kx_k^Hw$ ? Is there a formula ?

Q3: Equations : The final Equations are given by $w_{k+1}= w_k - \mu x_k({|y_k|}^2-1)y_k^T $ (T = transpose). While doing steepest ascent, we usually take a step in the opposite direction of the gradient. So, the minus should have become positive as according to the logic given in Least Mean Square.

I was looking at the Matlab implementation given here http://www.mathworks.com/matlabcentral/fileexchange/39482-blind-channel-equalization/content/lms.m
The Equations given in this impelemtation from lines 109 -- 113 are different from the ones in theory (excerpt above). The Equations in the implementation are :

for i=1:K
   e(i)=abs(c'*X(:,i))^2-R2;                  % initial error
   c=c-mu*2*e(i)*X(:,i)*X(:,i)'*c;     % update equalizer co-efficients
   c(EqD)=1;
end

What is correct? Can somebody please show the correct version?

As mentioned in the comments, the symbol $^H$ denotes the conjugate transpose of a matrix or vector, which means that the vector/matrix is transposed and that all of its elements are conjugated. The bar over $y_k$ means complex conjugate (note that $y_k$ is a scalar, not a vector or matrix).

When computing the gradient with respect to a complex variable in order to maximize or minimize a function, there's a trick which is explained in more detail in this answer. You basically take the derivative with respect to the conjugated variable and regard the non-conjugated variable as a constant. So when computing the gradient of the expression $\mathbf{w}^H\mathbf{x}\mathbf{x}^H\mathbf{w}$, you formally take the derivative w.r.t. $\mathbf{w}^H$, regarding $\mathbf{x}\mathbf{x}^H\mathbf{w}$ as constant, which yields the given result.

As for your third question, you should realize what you actually want to do. You want to minimize the error function, so you do not want to use a steepest ascent but a steepest descent algorithm, moving in the direction of the negative gradient.

And the Matlab code does in fact implement the equations you stated (replacing the modulus by the variable R2, which is equal to $1$ in your equations). Just use the error $e_k=|y_k|^2-1$ and the expression $y_k=\mathbf{w}^{(k)H}\mathbf{x}_k$ in your update equation, and you'll see it.


Difficulties acquiring sample period of a signal I'm teaching myself some DSP for future projects and my method is mostly learn by doing. I have this time-domain function with its correspondent frequency-domain transform:
$$s_1(t) = (t-2)e^{-t}u(t-2) $$
$$S_1(f)=\frac{e^{-2(1+j2\pi f)}}{(1+j2\pi f)^2} $$
Plotting the frequency-domain function in Matlab, for a frequency vector ranging from [-50,50) in steps of 0.1 (1000 values of frequency).

The question I don't know how to solve, is about the period sample. I know each sample is one of those dots, and that despite it can be arbitrary (choosed by us), there is criteria to follow in order to pick the best sampling rate possible to avoid aliasing.
If I set the frequency vector like I did, that means I already determined the sample rate, am I right? That also means, I can compute the sampling period as:
Ts=1/fs; fs=0.1 being the sample frequency.
Or should I go with another one? Should I try Nyquist (since its a continuous signal)? If so, how could I apply it?
Thanks for any suggestion you may have on this.

I think you are mixing up two things here. By "best sampling rate", I suppose your concern is how to select the $\Delta f$ in the spectrum such that the time-domain signal could be fully reconstructed from the sampled spectrum. This can be answered by the dual of the sampling theorem.

Consider the Fourier transform pair $s(t)\xrightarrow{\mathcal{F}}S(f).$ Let me first review the sampling theorem:

  Sampling in time with sampling rate $F_s=\frac{1}{\Delta t}$ causes replication in the frequency domain at every $F_s$ (i.e. frequency aliasing). To prevent frequency aliasing, $s(t)$ should be bandlimited to $B$, and the sampling interval should satisfy $$\Delta t=\frac{1}{F_s}<\frac{1}{2B}$$
  the signal's spectrum is recovered by filtering the replicated spectrum.

The dual theorem states that If the sampling interval in frequency domain is not short enough, it will result in time aliasing:

  Sampling in frequency with rate $\frac{1}{\Delta f}$ causes replication in the time domain at every $\frac{1}{\Delta f}$ (i.e. time aliasing). To prevent time aliasing, $s(t)$ should be time-limited to $T$ and the sampling interval should satisfy $$\Delta f<\frac{1}{2T}$$ The signal at the time domain can be recovered by time-windowing the replicated signal.

Getting back to your signal $s_1(t)$, since it is not time-limited (is not theoretically zero outside a time horizon), you should truncate it somewhere (e.g. at $T$). Then, it will satisfy the criteria and you can choose $\Delta f=\frac{1}{2T}$.


Coloring a waveform with spectral centroid or by other means Update 8

Unhappy of having to upload tracks to the service and looking at the new release of RekordBox 3 I've decided to take another look for an offline approach and a finer resolution :D

Sounds promising though still in a very alpha state:

Johnick - Good Time

Note that there's no logarithmic scale nor palette tuning, just a raw mapping from frequency to HSL.

The idea : now a waveform renderer has a color provider which is queried for a color for a specific position. The one you are seeing above gets the zero crossing rate for the 1024 samples next to that position.

Obviously there's still a lot to do before getting something robust but it seems like a good path...

From RekordBox 3:

Update 7

The final form that I'll adopt, much like in Update 3

(it's been a little Photoshopped to get smooth transitions between colors)

Conclusion is I was close months ago but did not consider that result thinking it was bad X)

Update 6

I've unearthed the project recently, so I thought about updating here :D

Song : Chic - Good Times 2001 (Stonebridge Club mix)

It's way better IMO, beats have a constant color etc ... it's not optimized though.

How ?

Still with http://developer.echonest.com/docs/v4/_static/AnalyzeDocumentation.pdf (page 6)

For each segment :

public static int GetSegmentColorFromTimbre(Segment[] segments, Segment segment)
{
    var timbres = segment.Timbre;

    var avgLoudness = timbres[0];
    var avgLoudnesses = segments.Select(s => s.Timbre[0]).ToArray();
    double avgLoudnessNormalized = Normalize(avgLoudness, avgLoudnesses);

    var brightness = timbres[1];
    var brightnesses = segments.Select(s => s.Timbre[1]).ToArray();
    double brightnessNormalized = Normalize(brightness, brightnesses);

    ColorHSL hsl = new ColorHSL(brightnessNormalized, 1.0d, avgLoudnessNormalized);
    var i = hsl.ToInt32();
    return i;
}

public static double Normalize(double value, double[] values)
{
    var min = values.Min();
    var max = values.Max();
    return (value - min) / (max - min);
}

Obviously there's a lot more code needed before you get here (upload to service, parse JSON etc) but that's not the point of this site so I'm just posting the relevant stuff for getting the above result.

So I'm using the first 2 functions of the analysis result, there is certainly more to do with but I still need to test. If I ever find something cooler than above I'll come back and update here.

As always, any hints on the topic are more than welcome !

Update 5

Some gradient using harmonic series

The color smoothing is ratio-sensitive otherwise it looks bad, will need some adjustment.

Update 4

Rewrote the coloring to happen at source and smoothed colors using an Alpha beta filter with 0.08 and 0.02 values.

A little better when zoomed out

Next step is to get a great color palette !

Update 3

Yellows represents mediums

Not so great when unzoomed, yet.

(the palette needs some serious work)

Update 2

Preliminary test using second 'timbre' coefficient hint from pichenettes

Update 1

A preliminary test using an analysis result from EchoNest service, note that it's not aligned very well (my fault) but it's much more coherent than the above approach.

For people interested in using this great API, start here : http://developer.echonest.com/docs/v4/track.html#profile

Also, don't get confused by these waveforms as they do represent 3 different songs.

Initial question

So far this is the result I get using a 256-samples FFT and calculating each chunk's spectral centroid.

Raw result of the calculations

Some smoothing applied (form looks much better with it)

Waveform produced

Ideally, this is how it should look like (taken from Serato DJ software)

Do you know what technique/algorithm I could use to be able to split the audio when average frequency changes over time ? (like the above image)

You can try first the following (no segmentation):

Process the signal in small chunks (say 10ms to 50ms in duration) - if necessary with a 50% overlap between them.
Compute the spectral centroid on each chunk.
Apply a non-linear function to the spectral centroid value to get a uniform distribution of the palette color used. Logarithm is a good start. Another option is to compute first the distribution of the centroid value other the entire file, and then color according to percentiles of this distribution (CDF). This adaptive approach guarantees that each color in the palette will be used in equal proportion. The drawback is that what is plotted in blue in one file won't be similar sounding to what is plotted in blue in another file if this adaptive approach is used!

It is not clear from the picture if Serato does this, or indeed goes one step further and attempts to segment the signal - it might not be surprising that it does since knowing the instants at which notes are present in a music signal can help synchronizing them! The steps would be:

Compute a Short-Term Fourier Transform (spectrogram) of the signal - using the FFT on short overlapping segments (start with a FFT size of 1024 or 2048 with 44.1kHz audio).
Compute an onset-detection function. I recommend you to look at this paper - the proposed approach is very effective and even a C++ implementation takes less than 10 lines. You can find an implementation in Yaafe - as ComplexDomainOnsetDetection.
Detect peaks in the onset-detection function to get the position of note onsets.
Compute the spectral centroid over all the time segments delimited by the detected onsets (don't forget to window and/or zero-pad!)
And don't forget the non-linear map! Note that the gradient-effect that appear between each note on the Serato waveform might be artificially generated.

Once you get this working, a "low hanging fruit" would be to compute a few other features on each segment (moments, a handful of MFCCs...), run k-means on these feature vectors, and decide on the color of the segment using the cluster index. See section II of Ravelli's paper.


Determining the resulting pixel value before creating a compressed JPEG image Is it possible to determine the resulting pixel values and modify them before the compressed jpeg image is actually created?

For example, I'm at the stage where I already have the quantized DCT coefficients of the original image and I want to change all pixels of the resulting image with a (255, 255, 0) value to (255, 254, 0).

What I was thinking is, if this is possible I would need to:

Get the equivalent DCT coefficient of original image and do quantization
Determine the resulting pixel values and do desired changes
Obtain the new DCT coefficients that will result to it
Proceed to Huffman encoding

Can the following be done during Huffman encoding?

Get the equivalent DCT coefficient from YCbCr and do quantizations
Do Huffman encoding and modify the desired pixel

To give you an upfront answer your question you want to determine the $pixel-values$ from $DCT-Coefficients$ probably there is nothing on earth better than doing IDCT itself!

However, leaving aside some vagueness i can re-phrase your question:

If i have an algorithm x - that is defined to work on Pixels, here is what i do:
DCT-Coeffs ->  IDCT  -> pixels  Algorithm_X  -> pixels  DCT  -> DCT-Coeffs
(ignoring other stuff that is peripheral).

  The question is: is there a short cut to this brute force where i am
  obviously undoing quite some stuff done in first stage?

What you are referring to is what is generally known as Compressed Domain Processing.

DCT-Coeffs ->  Compressed_Domain_Algorithm_X  -> DCT-Coeffs

Basically, all compressed domain algorithm exploits the linearity of DCT. For example

$$
Y = \alpha * X_1 + \beta * X_2
$$
then,
$$
DCT(Y) = \alpha * DCT(X_1) + \beta * DCT(X_2)
$$
where $\alpha$ $\beta$ are scalars.

There are many things possible, and while many NOT. It is too vast to answer the set of all compressed domain algorithms which are possible, but i would recommend key readings that will set a great stage for you.

Compressed Domain Processing of JPEG-encoded images by Brian C. Smith , Lawrence A. Rowe Real-Time Imaging (1996)
A Survey of Compressed Domain Processing Techniques by Brian C. Smith
Proc. NSF Workshop on Reconnecting Science
Algorithms for Manipulating Compressed Images, Brian C. Smith, Lawrence A. Rowe,
IEEE Computer Graphics and Applications, Sept. 1993, vol.13, (no.5):34-42.
Manipulation and compositing of MC-DCT compressed video by Shih-Fu Chang, David G. Messerschmitt IEEE Journal on Selected Areas in Communications, 1995/1, Volume
13 Issue 1
Compressed-Domain Video Processing by Susie Wee, Bo Shen, John Apostolopoulo
HP Labs Tech Report 2002.

Specifics:
If you are using libjpeg you can actually get quantized DCT co-efficients directly. However, most above literature refers un-quantized DCT co-efficients.


Musical onset detection approach I'm working on a real time beat detection algorithm. The first step, I believe, is to detect onsets (sudden volume increases). Then, in simplified terms, I can correlate these to find the tempo and use the strongest onsets to guess the beat position.

Say we are reading the audio in with 20ms chunks. My first guess was to just subtract the volume of one chunk from the previous one. The "volume" for a chunk would be a sum over the absolute value or squared signal, or maybe just the peak value (not sure which approach would be best). Then I can threshold these values to find the biggest volume changes.

To my surprise, according to papers like these (page 3) it seems that we often do onset detection by splitting into frequency ranges. For each audio chunk, they get the power for each of 7 frequency ranges, and compare it to the power of the previous chunk at that frequency (and the two neighboring frequency ranges), taking successive differences.

Then we have onset strengths, i.e. changes in power, for the 7 distinct ranges. And sum them back up across frequency to get the overall onset strength for that time.

Why is this done instead of just using the signal (no FFT required)? This seems to me to be just deconstruct and reconstruct the same information. The only thing I can think of is that the intermediate information of the power at different frequencies could be useful (isolating a kick drum, etc.)

Your question is really quite broad.

In a nutshell, if you think of the "successive differences", called spectral flux, understood per frequency range you can start to generalize to different kinds of music where maybe the kick drum is not present but the other instrumentation is sufficient to conclude that there is still on onset. One can also better measure the "strength" (loudness) of an onset, which is useful for algorithms processing onsets to conclude something at a higher rhythmic level (meter, phrasing, other groupings)

There's a lot of literature & code out there to try:

Look into MIREX publications and datasets
Try out code on github (https://github.com/CPJKU/SuperFlux, https://github.com/search?utf8=%E2%9C%93&q=onset+detection&type=Repositories&ref=searchresults, https://github.com/adamstark/Beat-Tracking-Evaluation-Toolbox)
Play with Beatseeker (https://www.ableton.com/en/packs/beatseeker/)


Term for signal waveform correction or demodulation Let say a periodic signal waveform $v(t)$ can be modeled as:

$v(t) = f(\phi) A \sin(2\pi ft + \psi)$

where $A$ is a constant amplitude, $f$ is signal frequency, $t$ is time and $\psi$ is a constant phase.

While $f(\phi)$ is $f = k \cos(\phi)$, where $k$ is a constant and $\phi$ is a function of time $t$.

Now the signal $v(t)$ is corrected by removing the dependency to $
\phi$ factor i.e. $v'(t) = v(t)/f(\phi)$.

My question: is the process getting $v'(t)$ from $v(t)$ considered signal waveform demodulation ? or undistortion ?

The generic term I can think safely would be signal waveform correction but it would be better if there is a precise terminology.

Demodulation is not accomplished by division...

So, your operation is undistortion, which has the better name signal restoration, in the signal processing context...


3D model from high speed videos? Is there a way to convert or somewhat replicate RADAR data from a set of given high speed videos taken from multiple angles.

The problem is quite ill defined but the task I have in hand is to reproduce RADAR data for a set of videos of soccer shootouts. I also have a corresponding file with the coordinates of the ball, speed of it but not quite sure how that will come into play with this.

You'd basically need to do two things:

find a 3D-model of your environment, and track the ball in that. This model needs to have knowledge of the materials involved. Since a video gives you no information on the involved object's radar reflectivity, scattering coefficients and other RF parameters, that's something that you'll have to enhance your model with – manually.
Simulate the environment (e.g. via ray-tracing/volumetric RF FEM) as RF environment, and in that simulation emit your radar signal and observe the echo.

Honestly, I don't think this is overly viable ask "task I have", unless "I" is very much more experienced in RF simulations than I am, or "I" is a medium to large team of researchers.


What to do for the right FFT in Matlab (two peaks and incorrect amplitude)? This is the FFT I'm using. But two peaks (did not expect the one in the beginning) are occurring and the amplitude is not as specified (expecting a value of 2). Any help is much appreciated!

Amp = 2;
freqHz = 10000;
fsHz = freqHz*2+1;
dt = 1/fsHz;

sine = Amp*sin(2*pi*freqHz*(0:dt:1-dt));

transform = fft(sine,fsHz)/fsHz;
magTransform = abs(transform);

faxis = linspace(0,fsHz/2,fsHz);
plot(faxis,fftshift(magTransform));
xlabel('Frequency')

The somehow symeetric aspect of the complex DFT might be misleading, so in the past I shared a simple Matlab function FFTR.m for real signals only, that might yield the result you'd expect ($2$ magnitude at peak frequency):

as long as you don't try bordeline endeavours. That requires a little normalization and a special care at DC and Nyquist. The code could be:

% Signal definition
Amp = 2;
freqHz = 10000;
fsHz = freqHz*2+1;
dt = 1/fsHz;
sine = Amp*sin(2*pi*freqHz*(0:dt:1-dt));

% Graphic part
[fftR,fftAxe] = FFTR(sine',dt);
plot(fftAxe,fftR,'o-');
xlabel('Frequency')
ylabel('Magnitude')


Determining the covariance of point clouds in real-time So basically I have a set of multidimensional data that I need to determine the covariance of between dimensions in real-time. Each point that comes in is a vector. I have gotten the mean and variance per dimension using this answer:
Determining the mean and standard deviation in real time. But I also need to find the covariance between dimensions, resulting in a covariance matrix. Each vector is 3 dimensional and is a point in space.

Does anyone know how to do this and does this?

To explain further, I want the covariance matrix so I can then use the eigenvalues and eigenvectors to characterize a point cloud. The point clouds are segments of a scanned environment using the Kinect. I need to make sure the segments are labelled consistently over frames. I will do this by matching the mean and covariance eigenvectors/eigenvalues of segments between frames. I am working with unity.

Edit:

Or what would be equally useful would be how to get the covariance from just the variances and means that I already have.

To make the Online variance algorithm
work for vectors, change one line:

M2 += np.outer( delta, delta2 )  # deltas 3-vecs, M2 3 x 3

# outer( x, y ) is pairs x[i] * y[j] --
# [[ x0 y0, x0 y1, x0 y2 ],
#  [ x1 y0, x1 y1, x1 y2 ],
#  [ x2 y0, x2 y1, x2 y2 ]]

A small python class for this is under gist.github.com/denis-bz .
It can covariance points with uniform weights (FIR),
or up-weight newer data (IIR), or both.

Is there a tutorial on your application with 3d point clouds ?


What kind of discrete convolution is this? Let x[n] and y[n] be defined for n=0,1, ..., N-1.

This (for example with Python's scipy fft / ifft)

a = ifft (fft(x) * fft(y))

should give the convolution

$$a[k] = (x * y)[k] = \sum_{n=0}^{N-1} x[n] y[k-n]$$

But is it done with

y[i] = 0 if $i<0$  ?

or

y[i] = y[N+i] if $i<0$ ?  (ex : y[-3] = y[N-3] )     (in this case do we call it circular convolution?)

Note : More generally, I have good books about signal processing / Fourier transform in general, but they are not really handy handbooks for discrete finite signals (x[n] for  n=0,1, ..., N-1). Do you have a good reference of handbook for this case (discrete finite signals) ?

Fast (FFT/IFFT) convolution results in circular convolution.

(But if sufficient zero-padding is explicitly added before the initial FFTs, circular and linear convolution would then produce identical results.)


Compressive Sensing Through MATLAB Codes I am new to the topic of compressed sensing. I read a few papers about it by R.Baranuik, Y.Eldar, Terence Tao etc. All these papers basically provide the mathematical details behind it, i.e., Sparsity, RIP, L1 norm minimisation etc. However, can anyone provide the MATLAB codes that achieve compressive sensing?

Thanks in advance for any help.

L1magic. This is the toolbox associated with the original paper.
CompSens. This looks like it's in C, but you could possibly call it with mex -- not sure.
Model-based compressive sensing toolbox.

  Most of the code is plain Matlab code
  Each folder in the package consists of a CS recovery algorithm based on a particular signal model, and a script that tests that recovery algorithm. The names of the scripts typically end with '_example.m'

Manopt. Possibly what you use to actually make the algorithms included in other toolboxes.

  Optimization on manifolds is a powerful paradigm to address nonlinear optimization problems.

Some examples from UMich.
Sparsity.

  This toolbox implements several algorithms to compute sparse expansion in redundant dictionaries and to solve inverse problems with sparse regularization (and also TV regularization).

But all of that, and more, is included in this list of toolboxes.

I've found the hard part is finding psuedocode -- that's where they actually describe the algorithm. Here are some examples of algorithms that included the psuedocode:

one paper that includes the actual algorithm (and the actual code available here in the function IST.
a beginners guide to compressed sensing.
Orthogonal matching pursuit
CS-KLT
Some more algorithms from Microsoft
MCMC for universal CS
An example
Compressive sensing of streams of pulses
GPU Accelerated Greedy Algorithms for Compressed Sensing
Adaptive sensing camera
Compressive sensing meets game theory


Low-pass filtering bounces when signal doesn't start at 0 Say I have an array representing a signal, going from 30 to -30 in 1000 increments, with a little bit of noise on those values. If I apply any king of low pass filtering (Besset, Butterworth, Savitsky-Golay, etc), then the 1st point is now very near zero and the next few points reach progressively for +30 after some wild oscillations depending on the filter and its parameters.

It's been a while (ahem, decades) since I last delved into the maths of DSP, so I'm a bit rusty and looking to know what I should do about it. I thought that the solution was to window my data before filtering it, but it doesn't seem to change much.

Obviously what I want to end up with is a nice 30 to -30 ramp, with the noise removed.

The problem is due to the fact that initially, the filter memory is filled with zeros (if not otherwise stated). Since the filter output is a weighted average of the filter memory and the input, it will start its output near zero, when the initial state is set to zero.

As Atul suggested, you can overcome this by padding 30 to the beginning of the sequence (to fill the memory before you actually filter the sequence of interest). However, the amount of padding depends on the filter and how fast it adjusts to the steady state (depending on filter order and e.g. Q-factor)

As a more elegant solution, you can initialize the filter state such that its first sample already outputs 30. You can calculate the filter state as the steady state response as to a unit step input. E.g. python has this functionality already built-in:

signal = np.linspace(30, -30, 1000)
sigma = 1
signal = signal + sigma * np.random.randn(len(signal))

plt.plot(signal)
b, a = scipy.signal.butter(5, 0.1)
zi = scipy.signal.lfilter_zi(b, a)

filtered = scipy.signal.lfilter(b, a, signal, zi=30*zi)[0]
plt.plot(filtered, lw=2)


price of dsp kit TMS320c6678(used for signal processing) This kit is used for implementing DSP concets in hardware
What is the meaning of highlighted price in attached fig?
Does it means 120 dollar per unit chip price?Figure

Yes. If you buy at least 1000 units, each will cost you 120


Impulse response from Frequency response: why using $e^{j\omega}$ as an input? Every resource that I can find uses this identity when deriving impulse response:

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

Suggesting that the input signal was $e^{j\omega}$. But by definition Impulse Response is the result of applying the filter to delta signal. Why then do we use $e^{j\omega}$ instead?

I think you misunderstand the notation. Writing $H(e^{j\omega})$ does not mean that $e^{j\omega}$ is the input signal. It just means that the frequency response is a function of the complex variable $e^{j\omega}$, because, as you might know, it is defined by

$$H(e^{j\omega})=\sum_{n=-\infty}^{\infty}h[n]e^{-jn\omega}=\sum_{n=-\infty}^{\infty}h[n]\big(e^{j\omega}\big)^{-n}\tag{1}$$

where $h[n]$ is the impulse response. Eq. $(1)$ is the discrete-time Fourier transform (DTFT) of the sequence $h[n]$. Generalizing to a complex argument $z$ you get the definition of the $\mathcal{Z}$-transform:

$$H(z)=\sum_{n=-\infty}^{\infty}h[n]z^{-n}\tag{2}$$

Comparing $(1)$ and $(2)$ shows that the DTFT is equal to the $\mathcal{Z}$-transform evaluated on the unit circle, i.e., for $z=e^{j\omega}$. Note that evaluating $(2)$ on the unit circle only makes sense if the unit circle is inside the region of convergence of $(2)$, otherwise $H(z)$ doesn't converge for $|z|=1$ and in that case the DTFT either doesn't exist or it takes a form which is different from $H(z)$ with $|z|=1$.

Note that some people just use $\omega$ as the argument of the function that expresses the DTFT, so you might as well come across the notation $H(\omega)$. After all, it's just a matter of convention.

Concerning the notation $H(e^{j\omega})$ you can also refer to this question and its answer.


Why is this MATLAB butterworth filter giving me wrong results Given a 128 Hz signal, in MATLAB I entered the following command

[B,A] = butter(2,[1/64,63/64],'bandpass'); %construct second order all pass filter
Filtered_Data = filter(B,A, Data) %filter data using butterworth filter

And got

This is far from what I expected since the filter should have passed all frequencies. an someone help point out what I'm doing incorrectly?

Your filter is passing only between 1 and 63 Hz, assuming that you meant to say that the sampling rate is 128Hz.

Since you are removing everything below 1Hz, the plot you have looks reasonable.

The range of values in your original signal is relatively small (4395 to 4420,) and the range of values in the output plot is also relatively small (cannot tell exactly due to scaling, but it looks small.)

The high values (+4000, -2000) in your output plot are the step response of your filter.  Your data jumps from zero (before the first data point comes in) to around 4400 at the beginning.  This step causes the filter to overreact - the short peaks up to 4000 and down to -2000 are the result.  This damps out right away, and the rest looks pretty much like what would be expected.

The offset (about 4400) looks like a DC component to the filter.  Since DC is lower than 1Hz, it gets removed as well.


How is compression ratio related to number of dct coefficients Lets say we have a picture with 6bits/pixel and we use discrete cosine transform with 8x8 blocks of pixels and each coefficient is 8 bits. We can use 2 8x8 masks, a zonal mask A and a zonal bit allocation mask B.

A = [1 1 1 1 1 0 0 0 1 1 1 1 0 0 0 0 1 1 1 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...... ]
B = [8 7 6 4 3 2 1 0 7 6 5 4 3 2 1 0 6 5 4 3 3 1 1 0 4 4 3 3 2 1 0 0 3 3 3 2 1 1 0 0 2 2 1 1 1 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0  ]

Is it correct to say that for A we have 64 : 15 compression namely 4.2 : 1 cause we keep 15 coefficients out of 64? does the same apply for B ?

Compression ratio is defined as the amount of bits in the original image divided by the amount of bits in the compressed image. In the original image the number of coefficients is 64, while you choose 15 coefficients after DCT in the compressed image with mask A. The ratio is 64*6/15*8. With mask B it is 64*6/41*8


Fourier components of $\cos(2\pi f_1t)$ I have the signal $s(t) = \cos(2\pi f_1t)$ and I am looking for its components vs the Fourier basis, over the interval $[0, T]$. The formula for computing the coefficients is
$$
s_n = \int_{t_0}^{t_1} s(t) \frac{e^{-j\frac{2\pi}{T} nt}}{\sqrt{T}} dt
$$
and in this case $t_0 = 0$ and $t_1 = T$. The solution seems to be
$$
\frac{\sqrt{T}}{2}(-1)^n\bigg[\mathrm{sinc}\big(f_1T+n\big)e^{-j\pi f_1T} + \mathrm{sinc}\big(f_1T-n\big)e^{j\pi f_1T}\bigg]
$$
but i can get no further than this
$$\frac{\sqrt{T}}{2}\bigg(\frac{e^{j2\pi (f_1T-n)}-1}{j2\pi (Tf_1-n)} + \frac{e^{-j2\pi (f_1T+n)}-1}{-j2\pi (Tf_1+n)}\bigg)$$

HINT:

Going from your last equation,
$$\frac{\sqrt{T}}{2}\bigg(\frac{e^{j2\pi (f_1T-n)}-1}{j2\pi (Tf_1-n)} + \frac{e^{-j2\pi (f_1T+n)}-1}{-j2\pi (Tf_1+n)}\bigg)$$
This can be simplified further down by considering the following:

\begin{align}
e^{j2\pi (f_1T-n)} &= e^{j\pi (f_1T-n)}\cdot e^{j\pi (f_1T-n)}\\
1 &=e^{j\pi (f_1T-n)}\cdot e^{-j\pi (f_1T-n)}\\
e^{-j\pi n} &= (-1)^n\\
\end{align}
And that
\begin{align}
\mathrm{sinc}(x) &= \frac{\sin(x)}{x}\\
\sin(x) &= \frac{e^{jx} - e^{-jx}}{2j}
\end{align}


digital butterworth filter design error I am trying to design a digital butterworth filter for the given specifications..
code is posted below..

rp=3;
rs=15;
FS=1;
wp=0.5*pi;
ws=0.75*pi;
pwp=2*FS*tan(wp/2);
pws=2*FS*tan(ws/2);
[n,wn]=buttord(pwp,pws,rp,rs,'s')
[b,a]=butter(n,wn,'s');
[bn,an]=bilinear(b,a,FS);%error
freqz(bn,an,512,FS);

My question is ,if i design a filter manually i am getting different values for 'bn' and 'an'..
but matlab gives different answer..
can anybody help me...?

I cannot run your code because I have not matlab on this pc, but I can try to give you some advice.

The first thing I will check is the way you defined the cut-off frequencies. The Matlab function uses normalized frequencies (look at the examples here http://it.mathworks.com/help/signal/ref/buttord.html), which means you have to define it in this way: lets assume you want to cut your signal at 5 Hz and your sampling frequency is $Fs$

cut_freq = 5;
Wn = cut_freq/(Fs/2);

Of course if you have a band-pass filter you have to define 2 cut-off frequencies and therefore $Wn$ will be a vector of 2 elements.

Hope this helps. As soon as I get my laptop I will try to run your code and see what happens.


extrapolate three points between line I have a line with known starting and ending 'x' and 'y' coordinates. I want to extrapolate three more between x and y points in MATLAB. I tried to find a lot but I couldn't get any help regarding extrapolate. Anyone please help me in this.

As you interpolate a line, I assume it is a linear interpolation. The solution reduces to a linear combination of the intial starting and ending coordinates. For instance, $i_1 = \frac{1}{3}s + \frac{2}{3}e $ and $i_2 = \frac{2}{3}s + \frac{1}{3}e $ provide the coordinates of two interpolation points $i_1$ and $i_2$ between start point $s$ and end point $e$.

clear
%% Initialization
pStart =[1 3]';
pEnd = [4 2]';
nPointInter = 3;

%% Interpolation weights (including pStart and pEnd)
lInterval = linspace(0,1,nPointInter+2);

%% Interpolation coordinates
pInterp = pStart*lInterval+pEnd*(1-lInterval);

%% Display
figure(1);clf
hold  on
plot(pInterp(1,:),pInterp(2,:),'xb')
plot([pStart(1) pEnd(1)]',[pStart(2) pEnd(2)]','or')
grid on

The result will be:


Impulse response of the inverse system to the backward difference I'm currently stuck on a problem. I have a DT LTI system $\mathcal{T}$ that depicts the first backward difference:

$$y[n] = \mathcal{T} \{ {x[n] \} } = x[n] - x[n-1]$$

Its impulse response is therefore: $h[n] = \delta[n] - \delta[n-1]$, where $\delta[n]$ is the DT unit impulse.

I want to compute the impulse response $h'$ of the system $\mathcal{T}'$ that, when cascaded with $\mathcal{T}$, recovers my input:

$$\mathcal{T}'\{\mathcal{T}\{x[n]\}\} = x[n]$$

I know the answer is that $h'[n]$ is the DT unit step $u[n]$, but I can't quite figure out why.

Here is my idea so far:

I know that the impulse response of the cascaded system is: $$h_{casc}[n] = \mathcal{T}'\{\mathcal{T}\{\delta[n]\}\} = \delta[n]$$
I know that this impulse response is the convolution of the impulse responses of $\mathcal{T}'$ and $\mathcal{T}$, so: $$\delta[n] = h[n] \ast h'[n] = h'[n] \ast (\delta[n] - \delta[n-1]) = h'[n] - h'[n-1]$$
So for $n=0$ I have: $$h'[0] - h'[-1] = 1$$
And for $n \ne 0$ I have: $$h'[n] - h'[n-1] = 0 \Leftrightarrow h'[n] = h'[n-1]$$

Still, this isn't enough to solve for $h'$ if I am not mistaken, as for example the following different signals would satisfy this:

$... 4 \ 4 \ 4 \ 4 \ 4 \ 5 \ 5 \ 5 \ 5 \ 5 \ldots$, Here: $h'[-1] = 4$, $h'[0] = 5$
$... 2 \ 2 \ 2 \ 2 \ 2 \ 3 \ 3 \ 3 \ 3 \ 3 \ldots$, Here: $h'[-1] = 2$, $h'[0] = 3$

Can anyone give me a hint what I am missing here or perhaps point me to an error or in the right direction?

You're on the right track. Let's call the impulse response of the second system (equalizing the first one) $h_2(n)$. The trick is to assume causality of the second system, i.e. $h_2(n)=0$ for $n<0$. By doing so, the solution for $h_2(n)$ becomes unique.

The total impulse response of the cascaded system is

$$h(n)=h_2(n)-h_2(n-1)\tag{1}$$

and this must equal $\delta(n)$. So we get

$$h_2(0)=1\\
h_2(n)=h_2(n-1),\quad n\ge 1$$

Which gives

$$h_2(1)=h_2(0)=1\\
h_2(2)=h_2(1)=1\\
h_2(3)=h_2(2)=1\\\vdots$$

So $h_2(n)=u(n)$, where $u(n)$ is the discrete-time step function.

Note that you can also find an anti-causal filter $h_2(n)$ satisfying the requirement $h(n)=\delta(n)$. This anti-causal solution is $h_2(n)=-u(-n-1)$. You can easily convince yourself that for this choice of $h_2(n)$ the total impulse response $h(n)$ given by (1) also equals $\delta(n)$.

Note also that all the solutions you gave at the end of the question have the form $h_2(n)=u(n)+c$ where $c$ is some constant. It is indeed true that all these systems are solutions to the problem because the frequency response of the first system has a zero at DC, compensating for any DC component of $h_2(n)$. However, all of these systems with $c\neq 0$ are neither causal nor stable.


Find Natural Frequency $\omega_n$ You are testing a new device to determine its characteristics. You find the magnitude of its frequency response shows a peak at 500 Hz, and the measured overshoot of the step response is 40% of the steady state level. Find the system’s characteristics and confirm your answers by generate plots of its frequency and step response using either Simulink or Pspice.

I'm trying to figure out how to find the natural frequency with the given info.

$\frac{500\,\text{Hz}}{140\,\%}\approx 357.145$;  I'm assuming my $r_o = 357.145$, where we settle at.

Given percentage overshoot (40%) I can solve for our damping ratio $\zeta = 0.279$.   I have formulas to solve for Natural Frequency ($\omega_n$) but not enough known variables to solve for it.  Does anyone know how to solve for $\omega_n$ here?

In our case, $\beta = 1.288$, if that helps.

Using this formula for the percentage overshoot, the damping ratio is given by

$$\zeta=\sqrt{\frac{\log(0.4)^2}{\pi^2+\log(0.4)^2}}=0.28$$

which is more or less the same value that you obtained. Consequently, you have an underdamped system. The peak frequency $\omega_{max}$ of the frequency response is related to the damping ratio and to the natural frequency in the following way (link):

$$\omega_{max}=\omega_n\sqrt{1-2\zeta^2}$$

With $\zeta=0.28$ and $f_{max}=500\text{ Hz}$ this results in a natural frequency $f_n=544.51\text{ Hz}$. Apart from a scaling factor, a second order system is completely characterized by $\zeta$ and $\omega_n=2\pi f_n$.


Is it possible to remove the non-random noise in the signal by averaging? Ex: We're measuring vibration Frequency response of a structure, and at the same time there is a constant source of vibration (noise) from a pump. The pump is exciting the structure by sinusoid correspond to its constant rotational speed.
Is it possible to remove the sinusoid which is induced by the pump using averaging technique (Signal averaging)?

Window-average filter is an example of lowpass filter. That means, if interesting frequencies in your signal are few times lower than noise, then yes, you can do that.
If they are close or higher, then no.
Use bandstop/highpass filter if your signal frequencies are much higher than noise, or adaptive filter like Recursive Least Squares instead, if your signal and noise frequencies overlap.
And don't forget, that your filter will affect your FRF measurements, so if you sample over the range that includes constant noise frequency, I'd suggest going with adaptive filter.
See also: What is the cut-off frequency of a moving average filter?


Zero-padding, energy conservation When calculating the FT of a signal, and the FT of that same signal but after having applied zero-padding to it, should't the amplitudes of the FT of the zero-padded signal decrease with respect to the FT amplitudes of the original signal?

1) Since the FT of the zero-padded signal has more frequency bins (since we have interpolated the spectrum), my intuition tells me that the amplitude of the frequency bins of this FT should be lower.

However, when I calculate the FT of the two signals (I use the fft function available in Scipy for Python), I get the same amplitudes for their spectra. I was not expecting this behaviour.

2) If I add the normalization term when calculating the FT (Scipy FFT does not include this), i.e. $FT = fft(signal) / length(signal)$, then I get different amplitudes for the two FT, which again, was not what I was expecting. Specifically, the zero-padded signal FT amplitudes decrease with respect to the non zero-padded, and the more zero-padding I add, the more the amplitudes of the FT decrease.

I am a bit confused with whether these results are supossed to be like that, or if I have done something wrong in my calculations.

What scipy.fftpack.fft effectively does is (docs):

y(j) = (x * exp(-2*pi*sqrt(-1)*j*np.arange(n)/n)).sum().`

This has no built-in normalization. If you do fft and then ifft, you need to normalize by multiplication by $1/n$ to get back your original data.

Or, you can do:

step 1: fft, then normalize by $1/\sqrt{n},$ then
step 2: ifft, then normalize by $1/\sqrt{n}.$

This gives the same final result because $1/\sqrt{n}\times1/\sqrt{n} = 1/n,$ and agrees with your intuition: The sum of squares of absolute values before and after each step remains equal. The steps compute the unitary discrete Fourier (DFT) and inverse discrete Fourier transforms (IDFT). With larger $n$ the same amount of "energy" as before will be distributed to more frequency bins so each value must be on average smaller in magnitude.


Subtracting Original Audio Signal here is my scenario: one speaker is sending a chirp, and a microphone is recording it as well as all following echoes. The function of the chirp is already know.

My goal: subtract the original chirp from the recorded sound, and leave only the echoes.

I think this questions is hard for me in these ways:

To match the original chirp and recorded sounds, I don't know which
approach to use: cross-correlation, or minimize sum(abs(difference
between two vectors)). The difference here means element-wise difference.
The echo may be mixed with the chirp itself, so I'd better only do
cross-correlation on the first small portion of the chirp.
The amplitude of recorded sound depends on the microphone.

Thank you very much!

Well my friend, let me answer your questions here. You've asked why subtraction of recorded sweep from original one will not produce any info about echos. So I took exponential sinusoid in range of $5\mathrm{Hz}$ up to $7995\mathrm{Hz}$, sampling frequency is $16\mathrm{kHz}$. Then I've used some crude filtering:

In result we get the following:

You can see that amplitude is modulated according to filtering applied. Now if we want to get some info about the impulse response of our system, then let's try simple subtraction as suggested by you, that's what we get:

What you can get from this time domain signal is what is the difference between different frequencies in time. Any info about reflections would not be straightforward. This is simple filtering, so there are none, but imagine you have some reverberation in the end, then your recorded signal is longer than played. Thus you have some tail of decaying sinusoids - nightmare to analyse anything
But when you deconvolve filtered signal with the inverse filter, that's what you get:

How beautiful is that! ;) You've obtained the impulse response of a filter that was applied at the very beginning. What's more, you can perform Fourier Transform of that impulse response and that will give you frequency response of this filter. I quickly did some frequency analysis, so you can see result below. Upper normalised frequency axis is forced to be logarithmic. You can easily relate that to filtering applied in the beginning. What's more, you can deduce that this plugin is using FIR filter because phase is linear!

If it was for example measurement of acoustic space, then your impulse response would be:

From that you can very easily detect any reflections, delays between arrivals via different paths and even flutter echo or different fancy stuff.

Because you are interested in detection of echo, then I suggest you to read the work of Dietsch (unfortunately in German): Ein objektives Kriterium zur Erfassung von   Echostörungen bei Musik- und Sprachdarbietungen

Above all following thesis: Evaluation of objective echo criteria. Guy is investigating echos a lot - you even have MATLAB files in the appendix.

Some bits about theory standing behind sweep-sine measurements: Advancements in Impulse Response Measurements by Sine Sweeps, and: Simultaneous Measurement of Impulse Response and Distortion with a Swept-Sine Technique

And bit more about derivatives from IR to convince you how powerfull it is: Theoretical and Applied Room Acoustics Including Parameters and IR Derivatives (although that one is targeting more at Architectural Acoustics)

I hope I did convinced you - good luck!


Discrete Fourier Transform: What is the DC Term really? I am currently toying around with the Discrete Fourier Transform (DFT) in Matlab to extract features from images. I like to fully understand the concepts that I use. I have read several explanations, such as this, but so far, none really explained the meaning of the "DC term". All I know is that the k'the term of the DFT can be written as:

where
 is the twiddle factor.

That means that the first term (the DC term), , is an amplitude without frequency.

Can someone explain why is it called the DC term? What is it's relation to "Direct Current"? And what is the relevance of the DC term? When is it useful, and for what?

The DC term is the 0 Hz term and is equivalent to the average of all the samples in the window (hence it's always purely real for a real signal). The terminology does indeed come from AC/DC electricity - all the non-zero bins correspond to non-zero frequencies, i.e. "AC components" in an electrical context, whereas the zero bin corresponds to a fixed value, the mean of the signal, or "DC component" in electrical terms.

As far as practical applications go, the DC or 0 Hz term is not particularly useful. In many cases it will be close to zero, as most signal processing applications will tend to filter out any DC component at the analogue level. In cases where you might be interested it can be calculated directly as an average in the usual way, without resorting to a DFT/FFT.


Clicks at end of chirp signal When creating a scipy chirp function I get some clicks at the start end the end.
The first one is quite easy to fix using a phase difference of 90 degrees.

The problem I have is that the end point is not 0.

for an static signal I can use the following code:

period_samples = fs / f
periods = np.ceil(T * fs / period_samples)
T = periods * period_samples / fs

and then I have the next zero point where I can stop.

But how can I do something similar with a chirp

my code now is:

import numpy as np
from scipy.signal import chirp

T = 1  # sec
fs = 44100  # samplerate
f0 = 20  # Hz
f1 = 20000  # Hz

t = np.arange(T * fs) / fs
Sig = chirp(t, f0, T, f1, 'logarithmic', 90)

From the code in Olli's comment, the phase is going to be

beta = t1 / log(f1 / f0)
phase = 2 * pi * beta * f0 * (pow(f1 / f0, t / t1) - 1.0)

So to ensure you have a zero crossing at the end of your signal, you need to choose a combination of t's end value and t1 for which the phase is modulo $2\pi$.

From your code, the last value of t will be 44099/44100. So the beta value is:

beta = 1 / log(1000)

and the phase of the last sample in the signal will be

phase =  2 * pi * 1/log(1000) * 20 * (pow(1000, 44099/44100) - 1.0)

which will be close to modulo $2\pi$ if you have a ratio of f1 to f0 equal to 1000.007748 instead of 1000 (after a bit of playing in Excel).  Try setting f1 to 20000.15495 instead of 20000.

EDIT : I made a mistake in the above and have corrected it. Here's a Google Doc that has the spreadsheet I've used (read only) for the correction.

Screenshot highlighting the close values below.


When I double bit rate I do not see a difference in I_spectrum and Q_spectrum, but why? I wanted to compare two different I-Q spectrum plots of a basic transmitter, one in which I doubled the bit rate. I expected a different result, that of possibly a higher/faster graph. It looks sort of the same, but that the doubled bit rate has more noise? I am not understanding the results here.
Bit rate of 1MHz

Bit rate of 2MHz

Do you know how ADS is plotting the spectrum? Plotting the spectrum without doing some kind of normalization will give you a higher magnitude. Once you determine the size of the FFT, normalizing by the length will give you the same magnitude no mater what sampling rate you choose.
For example let's take two rectangular signals, one sampled at 1 MHz and the other at 2 MHz. Below are their spectrums without normalization:

Since the bottom one is sampled twice as fast, it eventually produces an FFT size that is twice as long, hence the 6 dB increase in the peak.
Now compare this to the same exact signals, but now their magnitudes are normalized by their respective FFT sizes:

Now you can see that the peaks are the same magnitude. You can play with normalizing all day long to fit your need. It is the shape of the spectrum that is usually most important.
Here is some quick MATLAB code so you can maybe try it yourself and play around a bit.
%% Signal generation and FFT

% Sampling rates
fs1 = 1e6;
fs2 = 2e6;

% Rectangular pulse signals
t1 = 0:1/fs1:1e-5;
t2 = 0:1/fs2:1e-5;
pulseSignal1 = ones(1, numel(t1));
pulseSignal2 = ones(1, numel(t2));

% FFT setup
nfft1 = 100*numel(t1);
f1 = fs1.*(-nfft1/2:nfft1/2-1)/nfft1;

nfft2 = 100*numel(t2);
f2 = fs2.*(-nfft2/2:nfft2/2-1)/nfft2;

%% Without Normalization

figure;
subplot(2, 1, 1);
plot(f1./1e6, 20*log10(abs(fftshift(fft(pulseSignal1, nfft1)))));
xlabel("Frequency (MHz");
ylabel("Magnituide (dB)");
legend("F_s = 1 MHz");
ylim([-40 50]);

subplot(2, 1, 2);
plot(f2./1e6, 20*log10(abs(fftshift(fft(pulseSignal2, nfft2)))));
xlabel("Frequency (MHz");
ylabel("Magnituide (dB)");
legend("F_s = 2 MHz");
ylim([-40 50]);

%% With Normalization

figure;
subplot(2, 1, 1);
plot(f1./1e6, 20*log10(abs(fftshift(fft(pulseSignal1, nfft1)./nfft1))));
xlabel("Frequency (MHz");
ylabel("Magnituide (dB)");
legend("F_s = 1 MHz");
ylim([-80 -10]);

subplot(2, 1, 2);
plot(f2./1e6, 20*log10(abs(fftshift(fft(pulseSignal2, nfft2)./nfft2))));
xlabel("Frequency (MHz");
ylabel("Magnituide (dB)");
legend("F_s = 2 MHz");
ylim([-80 -10]);


Equalizers with low data rate systems In case of using equalizers with ZigBee systems, what type of fading does the channel show, flat or frequency selective, fast or slow?
Is the inter-symbol interference (ISI) problem encountered in ZigBee systems?

I am asking this question because I konw that ZigBee Systems transmit data at low rate (250 kb/sec), and ISI problem and hence frequency selective fading appear with high data rate applications. I found some papers employing equalizers with ZigBee under Rayleigh fading conditions, but the type of fading was not indicated.

In the $2.4 \textrm{ GHz}$ band, the data rate is indeed just $250\textrm{ kbit/s}$ but the direct-sequence spread spectrum (DSSS) coding results in a bandwidth of $2\textrm{ MHz}$. This means that there will be flat fading only if the channel's coherence bandwidth is significantly greater than $2\textrm{ MHz}$, say around $10\textrm{ MHz}$ (this is just a fuzzy estimate). This corresponds to a maximum delay spread of approximately $1/(10\textrm{ MHz})=100\textrm{ ns}$.

Consequently, any environment with a delay spread greater than approximately $100\textrm{ ns}$ will result in frequency-selective fading. So in relatively small indoor areas the fading will be mostly flat, but in large office spaces or outdoors you can expect frequency-selective fading.


Is Kelly Lochbaum Ladder Filter LTI? I am a bit confused about the Kelly Lochbaum Ladder Filter. Image below taken from here.

My question is, is this an LTI (linear time-invariant) filter? I thought that it was LTI, but I am not sure anymore.

As long as the graph consists only of delays, add/subtract and multiply-with-a-constant, then the system will be LTI.

In this case, if the coefficients $k_1 ... k_m$ don't change with time, the system is LTI.


Zero State and Zero Input Responses from Steady State Response Let a linear, time-invariant, causal discrete linear system be described by the following difference equation:

$$y[n] = 0.9y[n-1] + 0.1x[n]$$

Assuming that $y[-1]=2$ and $x[n]=20\cos(\Omega n)u[n]$, $\Omega = \omega T_s = 0.2\pi$, find the total response, identifiying the terms related to the zero state and zero input responses.

I've found the transfer function of the system to be:
$$H(z) = 0.1\frac{z}{z-0.9}$$

and from it I've determined the magnitude and phase of the steady-state response:
$$|H \left(e^{j \Omega}\right)| = \frac{0.1}{\sqrt{1.81-1.8\cos(\Omega)}}$$
$$\angle H \left(e^{j \Omega}\right) = -\arctan\left[\frac{-\sin(\Omega)}{1-0.9\cos(\Omega)} \right]$$

Therefore the steady-state response of the system is:
$$y_{ss}[n] = 20 \frac{0.1}{\sqrt{1.81-1.8\cos(\Omega)}} \cos\left(\Omega n -\arctan\left[\frac{-\sin(\Omega)}{1-0.9\cos(\Omega)} \right] \right)$$

$$\therefore \boxed{y_{ss}[n] = 3.363\cos(0.2\pi n + 1.138)}$$

But how do I find the total response (zero state and zero input) from the steady-state response?

Some books prefer different notations, but I will assume that zero-state response is the one due to pure input $x[n] = \cos(0.2 \pi n) u[n]$ under system zero-initial conditions, while the zero-input response is the one with $x[n] = 0$ and initial state being $y[-1] = 2$.

Now, first the zero-input ($x[n] = 0$ ) response, purely due to the nonzero initial condition ($y[-1] = 2$), of the system can be found by a simple recursion of the causal system as:

$$
\begin{align}
y[0] &= 0.9 \cdot 2 \\
y[1] &= 0.9^2 \cdot 2 \\
y[2] &= 0.9^3 \cdot 2 \\
... &= ... \\
y[n] &= 2 \cdot 0.9^{n+1} u[n] \\
\end{align}
$$

which yields for the zero input response (for $n \geq 0$) as

$$\boxed{ y_{zi}[n] = 2 \cdot (0.9)^{n+1} u[n] }$$

Now we shall find the zero-state response, assuming zero initial conditions. This will be lengthy because of the input form $x[n] = 20 \cos(0.2 \pi n) u[n]$, we can follow the following approach.

Let $X(w)$ be the Fourier transform of the input, and $H(w)$ be the Freuqency response of the system described by the above differential equation, then the output $y[n]$ will have the Fourier transform given by $$Y(\omega) = H(\omega) X(\omega).$$

Now $$H(\omega) = \frac{0.1}{1- 0.9 e^{-j \omega} } $$ and you have to find out what $X(w)$ is, then multiply them to find out $Y(\omega)$ and perform an inverse Fourier transform, after a simple partial fraction expansion. Denoting this solution as $y_{zs}[n]$, then your total solution is

$$y[n] = y_{zi}[n] + y_{zs}[n]$$

Now, try this on your own and see if you can complete the solution... Then write a matlab program to check your result.


Alignment of 2 Set of Samples from Different Sensors If we measure the heart rate of a subject with two different devices, which have big different sampling rates then how we can compare their outcomes. For instance, one of the devices has the sampling rate of 1000 while the other one has sampling rate of 32.

There are 2 main factors here:

Aligned timeline
Namely each sample has time stamp taken with the same reference time.
The lower sampling rate is high enough according to the sampling theorem
Namely even the lowest sampling rate, in your case 32 [Hz] is high enough to reconstruct the original signal.

If both (1) and (2) hold you'll be able easily to compare them by interpolating both into a shared time grid.
If (1) doesn't hold you may be able, by interpolation and some other DSP techniques (Such as Matched Filter and Dynamic Time Warping), to align both to the same grid.
Yet it would be harder (And in some cases, not achievable).


Transfer function given as z-Domain fraction: gain calculation I am stuck on a question regarding transfer functions. I have the answer and have attempted the question but am having trouble solving it.

Q: Given the transfer function of a digital filter is (5z-5)/z, compute the gain (in dB) of the digital filter at the cutoff frequency (10 rad/s)

The answer given is 14.65dB. (or a gain of 5.401)

I am unsure how to do this:

I see that the transfer function can be rewritten as 5(z-1)/z leading me to believe the gain would be 5 which is 13.98 dB, more than half a dB off.

What am I doing wrong here?

Obviously, you might need to use the given frequency somehow to solve the problem. Guessing alone will usually not suffice. Since this is a homework type problem, and since you still need to learn a lot of basics, I'll give you a few hints to get you started.

First, the frequency domain behavior can be obtained from a (stable) transfer function by substituting $z=e^{j\omega T}$, where $T$ is the sampling period, which equals the inverse of the sampling frequency $f_s$ ($T=1/f_s$). So you won't be able to solve the problem without knowing the sampling frequency (unless the frequency is given as a normalized value, which is not the case in your exercise).

The gain at the given frequency $\omega_c$ is then given by the magnitude of the transfer function evaluated at $z=e^{j\omega_cT}$:

$$g=\big|H(e^{j\omega_cT})\big|$$

I trust that you know how to convert that linear gain to the corresponding value in dB.


Insertion loss equality - proof From this paper (slide 11), the insertion loss (or attenuation) is defined as

$$
L(\omega^2)=\frac{\lvert V_i(j\omega)\rvert^2}{\lvert V_o(j\omega)\rvert^2}=\frac{1}{\lvert H(j\omega)\rvert^2}=10 \log\left(\frac{1}{H(j\omega)H(-j\omega)}\right)
$$

I'm interested specifically in, where does the last inequality follow?

As pointed out in PeterK.'s comment, the way it is written this does not make sense. What is meant here is that the insertion loss can be measured on a linear scale (i.e., as a factor) as

$$L(\omega^2)=\frac{1}{|H(j\omega)|^2}$$

or it can be given in decibels:

$$L_{dB}(\omega^2))=10\log_{10}(L(\omega^2))$$

The remaining question could be why $|H(j\omega)|^2$ can be replaced by $H(j\omega)H(-j\omega)$. The reason is the symmetry of the frequency response $H(j\omega)$ for real-valued filters:

$$H(j\omega)=H^*(-j\omega)$$

From this it follows that

$$|H(j\omega)|^2=H(j\omega)H^*(j\omega)=H(j\omega)H(-j\omega)$$


Minimum Mean Square Estimator - Equivalent Expressions to Minimize Given $ M \in \mathbb{R}^{N \times N} $ which is a Positive Definite Matrix.
Let $ \hat{x} $ the MMSE of $ x $ given $ z $, namely $ \hat{x} = \mathbb{E} \left[ x \mid z \right] $.

Prove the equivenalce of the following expressions:

$ \arg \min_{\hat{x}} \mathbb{E} \left[ {\left( \hat{x} - x \right)}^{T} \left( \hat{x} - x \right) \mid z \right] $.
$ \arg \min_{\hat{x}} \mathbb{E} \left[ {\left( \hat{x} - x \right)}^{T} M \left( \hat{x} - x \right) \mid z \right] $.
$ \arg \min_{\hat{x}} \operatorname{Tr} \left( M \mathbb{E} \left[ {\left( \hat{x} - x \right)}^{T} \left( \hat{x} - x \right) \mid z \right] \right) $

Equivalence means all are actually minimized by $ \hat{x} = \mathbb{E} \left[ x \mid z \right] $.

The equivalence of 2 and 3 is easy using the Cyclic Property of the $ \operatorname{Tr} \left( \cdot \right) $ operator.
Yet showing the invertible transformation keeps the solution (Equivalence of 1 and 2) isn't trivial.

Since $ M \in \mathbb{S}^{N}_{++} $ (In other convention $ M \succ 0 $) by Cholesky Decomposition there is a Triangular Matrix $ R \in \mathbb{R}^{N \times N} $ such that $ M = {R}^{T} R $.

Using this fact one could prove $ 1 \iff 2 $ as following:

$$\begin{align*}
\arg \min_{\hat{x}} \mathbb{E} \left[ {\left( \hat{x} - x \right)}^{T} M \left( \hat{x} - x \right) \mid z \right] & = \arg \min_{\hat{x}} \mathbb{E} \left[ {\left( \hat{x} - x \right)}^{T} {R}^{T} R \left( \hat{x} - x \right) \mid z \right] && \text{$M = {R}^{T} R $} \\
& = \arg \min_{\hat{x}} \mathbb{E} \left[ {\left( R \left( \hat{x} - x \right) \right)}^{T} \left( R \left( \hat{x} - x \right) \right) \mid z \right] && \text{} \\
& = \arg \min_{\hat{x}} \mathbb{E} \left[ {\left( R \hat{x} - y \right)}^{T} \left( R \hat{x} - y  \right) \mid z \right] && \text{Defining $ y = R x $}
\end{align*}$$

Clearly, the above, as the classic MMSE problem, with respect to $ y $ given $ z $ which is minimized when:

$$ R \hat{x} = \mathbb{E} \left[ y \mid z \right] \Rightarrow \hat{x} = {R}^{-1} \mathbb{E} \left[ y \mid z \right] $$

Since $ M \succ 0 $ then $ {R}^{-1} $ is defined and the above is valid.
Moreover by linearity of the Expectation Operator:

$$ \hat{x} = {R}^{-1} \mathbb{E} \left[ y \mid z \right] = {R}^{-1} \mathbb{E} \left[ R x \mid z \right] = {R}^{-1} R \mathbb{E} \left[ x \mid z \right] = \mathbb{E} \left[ x \mid z \right] $$

Hence the minimizer of 1 indeed minimizes 2.

Using the cyclic property of the Trace Operator and the linearity of the Expectation Operator (Hence one could change the order of the $ \operatorname{Tr} \left( \cdot \right) $ and $ \mathbb{E} \left[ \cdot \right] $) one could easily show $ 2 \iff 3 $ hence the problem is solved.


Does only filter coefficients can represent it's own filter characteristic? I'm wondering that does only filter coefficients can represent it's own filter characteristic?

I found filter design tool as FILTERDESIGNER of matlab. that can make filter coefficients of low pass filter or band pass filter or high pass filter etc as the below.

But I want to know what if there is only information of coefficients without information which is Response Type such as low or high and band filter.

Can we know that coefficients was made originally for low pass filter with  h[k]*x[n-k] equation?

You can compute the filter's frequency response to check the behavior of the filter:

$$H(\omega)=\sum_{n=0}^{N-1}h[n]e^{-jn\omega}\tag{1}$$

where $h[n]$ are the filter coefficients, and $N$ is the filter length (i.e., the number of coefficients). The frequency response $H(\omega)$ is a complex function. You can compute its magnitude and phase to evaluate the filter's magnitude and phase characteristics.

You can compute (and plot) the frequency response using Matlab's freqz command.


Why we measure bandwith at -3dB? To characterize a signal it's a common use to determine bandwidth at -3dB from the peak (half power). Also, filter's cutoff frequency is pointed at -3dB. Why we use this rule?

It is a common value because at $-3 \ \mathrm{dB}$ the power of the signal is reduced to half its value. I'll show a brief example to make it clearer.

Suppose that you have a signal whose amplitude is $1 \ \mathrm{V}$. If we reduce its amplitude by $-3 \ \mathrm{dB}$, the new amplitude would be close to $ 0.707 \ \mathrm{V} $.

The power of a signal is proportional to the square of its amplitude. For simplicity, let's assume a resistance of $1 \ \Omega$ on which we are going to calculate the dissipated power. Therefore, originally we would have had a power of $$(1 \ \mathrm{V})^2\cdot\frac{1}{1 \ \Omega}=1 \ \mathrm{W}$$

But after the reduction, we have a power of

$$(0.707 \ \mathrm{V})^2\cdot\frac{1}{1 \ \Omega}=0.5 \ \mathrm{W}$$

The power has decreased to half its original value. That is why the $-3 \ \mathrm{dB}$ point is commonly used to characterize frequency-related properties.


Processing and plotting quadrature data in the frequency domain Ok, so I recently started playing around with an SDR in an effort to challenge and educate myself. So, let's say the sampling rate (Fs) is 3MHz and of those we grab 3M samples. Now we are left with 3M samples consisting each of an IQ pair. Complex data. My understanding is that I would perform an FFT using the entire sample, consisting of IQ as the input for the FFT which would yield another complex answer. From this, calculate the magnitude and plot starting from 0? I'm not entirely positive which samples to throw out, which ones to keep, and in which order to plot them.

Do I perhaps separate the real and imaginary prior to the FFT for the first 1.5M samples, perform an FFT for 1.5M of I and FFT for 1.5M of Q and plot the real to the right of 0 and the imaginary to the left of 0? Throwing out the rest of samples from 1.5-3M?

Any help is very much appreciated. Thanks all.

You do not need to break it into real and imaginary. If you are interested in the spectrum, you will take the FFT and look at its squared absolute value i.e., magnitude square. Better see it in log domain, i.e., $20\log_{10}(|X(f)|)$, starting from 0 to number of FFT points.

For 1Hz resolution you need to take FFT of all 3M samples and FFT size should be 3M. But to ease the computational load you can consider taking a window of 1024 samples then take its  FFT using an FFT size of 1024 (In that case each FFT bin corresponds to 3M/1024 Hz). The number 1024 is just an example, you can consider any other value $N$ that suits you.

For plotting take the magnitude response of each window and plot it. This is spectrum plot without any overlap.

For 50% overlap move the window by $N/2$ samples, take another mag. resp. It is better to average the response over few windows.

Windowing is another interesting parameter to play with. In the above case the window is rect, i.e., $w(k) = \sum_{n=0}^{N-1} \delta(k-n)$, where size of window is N. If input signal is $x(k)$ then to pick N points we simply multiply $w(k)$ by $x(k)$. You can see other windows like Kaiser, Hamming or Hanning, and use which ever suits your needs.


Is subtraction a common method for efficient filtering? I've recently been doing some DSP programming with regard to filter banks. One topic that I've come across was that of the dyadic filter bank. This filter bank essentially breaks a signal into sub-band signals corresponding to some octave.

Looking at some diagrams, I notice that many of them illustrate implementing both a highpass filter and low pass filter. This will essentially generate two sub-band signals. Conceptually, this makes sense.

However, it seems rather inefficient to implement two separate filters. Would it not be more efficient to implement a single filter, and then subtract the output from the original signal to determine the other signal?

In other words, implement just the low pass filter to obtain the lower subband signal. Then simply subtract that lower subband signal from the inputted signal. This should yield the higher band signal.

So my question:

Is this technique common?
Furthermore, would it be vulnerable to phase distortion in the sense that your low pass filter may have phase distorted or lagged the lower band signal?

You are right about the importance of the phase.

Subtraction in the time domain is equivalent to subtraction in the frequency domain. If the frequency domain phase at some frequency is not equal between the complex number being subtracted and the complex number subtracting, even if their magnitudes are the same, the result will be non-zero:

Figure 1. Subtraction of complex numbers of equal magnitude does not always result in a zero. The calculation is black = red - blue.

If the black dot represents the combined frequency response for a system output = input - filtered input at a frequency that you intend to be part of a stop band, you must ensure that the filter is zero-phase. Or if the filter is linear-phase and has in its impulse response a center of symmetry at some time, you must delay the copy of the input by the same amount of time before subtraction. This will not work with efficient infinite impulse response (IIR) filters because of their non-linear phase.


suitable application for 2d FFT I have been building a convolution library for multichannel audio and am using FFTW as the main FFT transform library. There are 4 channels of audio in play at any one time and they are interleaved, so one "frame" contains 4 separate samples - one from each channel.

I have been using the Real to Complex 2D functions up until now, but I am starting to be concerned that multichannel audio is perhaps not the right application for this - even though it fits with the documented data structure.

Should I be using 4 separate 1D transforms instead of a single 2D FFT pass?

You should consider what you're doing.

Let me draw your audio data in the manner I understood it's in your memory. I'll call your channels A,B,C,D, and number their respective samples increasingly:

A0 B0 C0 D0 A1 B1 C1 D1 A2 B2 C2 D2 …

Now, if you want to apply a 2D FFT, your data needs to have the structure of 2D data – i.e. a matrix. You didn't quite specify, but I assume you've got a matrix

$$
\mathbf S = \begin{pmatrix}
A_0 & B_0 & C_0 & D_0\\
A_1 & B_1 & C_1 & D_1\\
A_2 & B_2 & C_2 & D_2\\
A_3 & B_3 & C_3 & D_3\\
 \vdots & \vdots & \vdots & \vdots\\
A_{N-2} & B_{N-2} & C_{N-2} & D_{N-2}\\
A_{N-1} & B_{N-1} & C_{N-1} & D_{N-1}\\
\end{pmatrix}
$$
And you're applying the $4\times N$ DFT (FFT):

$$\begin{align}
\hat{\mathbf{S}} &= \text{DFT}_{4\times N}\left\{\mathbf S\right\}\\
&= \left(\frac1{\sqrt{4N}}\sum_{n=0}^N\sum_{m=0}^4 S_{n,m}\, e^{j2\pi\left(\frac{nk}N + \frac{ml}4\right)} \right)_{k,l}\tag1\\
&=\left(\frac1{2\sqrt N}\sum_{n=0}^N
A_n e^{j2\pi\left(\frac{nk}N + 0\right)}
+ B_n e^{j2\pi\left(\frac{nk}N + \frac l4\right)}
+ C_n e^{j2\pi\left(\frac{nk}N + \frac {2l}4\right)}
+ D_n e^{j2\pi\left(\frac{nk}N + \frac {3l}4\right)}
 \right)_{k,l}\\
&=\left(\frac1{2\sqrt N}\sum_{n=0}^N
A_n e^{j2\pi\frac{nk}N}
+ B_n e^{j2\pi\frac{nk}N} e^{j\frac12\pi l}
+ C_n e^{j2\pi\frac{nk}N} e^{j\pi l}
+ D_n e^{j2\pi\frac{nk}N} e^{j\pi \frac32 l}
 \right)_{k,l}
\end{align}
$$

Looking at the last line: you've got to ask yourself whether you really want to have the DFT over each of your 4-channel frames – that sounds intuitively wrong.

Look at eq. $(1)$: you can swap the summations; hence, it doesn't matter if you first take the FFT of every channel column vector in "time direction", and then the row-wise FFT of the resulting matrix in "channels direction", or first do the row-wise FFT over all four channels, and then over all times.

I think you really want 4 1D-FFTs of length $N$, and not one $N\times 4$ 2D-FFT.

Luckily, the advanced 1D-FFT real-to-complex API of FFTw makes it really easy to do four identical FFTs on interleaved data. I annotated the API with the values that seem to be apporpriate for your use case:

fftw_plan fftw_plan_many_dft_r2c(int rank, // Dimensionality == 1
        const int *n,       // FFT input length == N
        int howmany,        // Amount of FFTs == 4
        double *in,     // pointer to A0
        const int *inembed, // pad/truncate input ? No --> NULL
        int istride,        // Distance between input elements that belong to the *same* FFT == 4
        int idist,      // Distance between the first elements of different FFTs == 1
        fftw_complex *out,  // pointer to output space
        const int *onembed, // see above inembed
        int ostride,        // same as above
        int odist,      // same as above
        unsigned flags);    // allow destruction of input buffer? Measure or guess the best algorithm?