Malaysia Flight 370 Image Processing I wish that all those people on that flight were with us, and maybe they still are, as I like to think they are on an island somewhere waiting to us to find them ... however, it does not look good...

Recently, a spy plane searching for possible remnants of Malaysian Flight 370 took this photo:

However they cannot tell what this is an image of, floating in the sea.

I was wondering, what image processing techniques might us in the DSP community utilize here to better ascertain / clean up the image so that we might be able to get a better picture?

The waves of the ocean would seem to contribute to high-frequency components, so perhaps nulling the high frequency fourier co-efficients in the 2D-FFT, and then inverse transforming?

What other techniques might we use? ... perhaps we can segment it while respecting shapes of big features?...

There's little chance of getting any useful info from that one photo.

What would help is to have a series of photos.  With a series, even all from the same angle, one could average out the moving waves and align the suspected object.  It would be interesting to see if the object is tilting, bobbing up and down, or otherwise moving.

Images from different angles would help determine the 3D shape of the object, and determine if it's on the water surface, partially submerged, or entirely under but near the surface.

There's unlikely to be any interesting image processing work here, just routine contrast adjustment, some smoothing out of noise, maybe a bit of sharpening.


modelling wireless multi-path channel - is channel response sum of delta(t-path_delay) or "delta" should be smoothed somehow? Consider the wireless channel. It is mult-path - means that signal is sum of sent signals with some delay and fadings.
The standard models for "impulse response" for multi-path channel in wireless is:

$ h(t) = c_1 \delta(t - t_1) + c_2 \delta(t - t_2) + ... + c_N \delta(t - t_N) $,

where $t_i$ are path-delays and $c_i$ are fadings.

Question I wonder is this model realistic  ? Is it really reasonable to assume "delta"-functions for each path or may be it we need to assume something more smooth ?
What is good reference, where the results of real measurements can be found ?

[EDIT] in response to comment.

I am interested in cellur telephone links. I know that the standard e.g. channel 3GPP models assume 6-9 paths with delta-functions. I am  suspicious about this assumption. So I want to see results of real measurements which support this model.

[EDIT 2] More details.
I am interested in channel estimation for OFDM system cellur system. Assume carrier frequency is something like 1GHz - 2GHz.
One user typically may have ~100 OFDM elements. Each  element is about 15kHz.
So the user bandwidth 1.5MHz.

What am suspicious about is that if we have 9 paths and we have about 100 pilot symbols, then  our channel estimation would be quite good , and BER will be close to ideal ch.est. Because we need to estimate ~9*2 unknowns from 100 of data - so we get good quality of estimation. This is in theory, but I am afraid if will test such system in practice
we mail fail, cause in practice it might be not 9 paths but 90 paths ...

So let me stress the following part of the question:

Where to find results of real measurements in urban requirement for  channel response ?

I guess that should be old stuff may be from 60-ies ? Cause many modern literature concern MIMO channels and pay attention to angular spread and so on - while I need more simple information.

Depends on the physics of the situation. The model is reasonably good if you have secular reflections and it is not very good if you have diffuse reflections. If the size and geometry details of the reflectors are large as compared to the wavelength, than its mostly specular. If the surfaces are uneven are very non flat, then its more diffuse. This is all measured in units of wavelengths. For a 5GHz link that's about 6cm. For an FM radio link that's about 3m.
It also depends on the relative bandwidth of the signal. For most wireless links that's pretty narrow, i.e the behavior at the carrier frequency is a good model for the behavior of the overall system. Most natural physical transfer functions such as a reflection of a wall do not change very rapidly with frequency so that's why the singe impulse is a good starting point.
Of course there are always exceptions to the rules so it depends on your specific problem.


Sampling the Wave Equation for String Vibration I am having a difficult time understanding the basic principles when sampling the function describing the vibration of a string. Mainly I am confused on given a sampling rate how many samples to take and what size the increments when sampling should be.
I am using the solution to the wave equation from this article (in german).
If my explanation is unclear let me know and I'll try to give more information.

Assuming I understand the question correctly, you need to sample the string's position at least twice as fast as the maximum frequency you anticipate the string can move or the highest frequency you care about. So if the highest overtones the string produces are 16kHz, you need to measure at at least 32kHz if you want to be able to reconstruct the original motion faithfully. This is the Nyquist-Shannon sampling theorem and you can probably find several hundred descriptions of it here and elsewhere on the internet.
If you know the lowest frequency the string can produce or that you care about, there are other tricks you can do (look up heterodyne) to reduce your sample rate requirements assuming you have the equipment to do so on the analog side.


Beat detection in Matlab or Python Try to determine the start time and duration of each "beat" in an audio signal (circled in figure below). It looks like a simple question, and python or Matlab should have toolbox/functions to do this, just couldn't find them. I tried low pass filtering and freq analysis, but neither worked.

The area you're asking about is called "Onset Detection" this is used all throughout signal processing, but most oven used in Music Information Retrieval to identify the start and stop to notes.

While there are many ways to accomplish this lets consider two: energy and  frequency domain.

Energy: We know that when there is a spike in a signal there is a marked increase in the energy of the signal. You could simply create sampling windows and calculate the average absolute energy of the window by tracking the area under the curve. When you see a large increase it means you've reached a peak (or close enough to an onset)

Frequency:
We know that in the time domain, if a signal is high in amplitude and short in time it will have a wide spread in the frequency domain with most of that content in the higher end of the frequency spectrum. (Dirac function is a prime example of this)

You could cut your signal into small windows, calculate the FFT and detect when there is a large change in high frequency content.

These two methods are theoretically the same, however your implementation would result in different accuracy.

Here's a tutorial on many of the onset detection algorithms out there.


Discreteness and periodicity in Fourier transform Why discreteness in time / frequency domain dictates periodicity in the other frequency / time domian?

For example the DTFT is perodic in frequency?
Why it doesn't contain all the frequencies?
Why the DTFS contain only finite number of frequencies?

NOTE:

I've read the derivation but sometimes i see it's like a definition or a claim in the case of DTFS.

How we prove that equation of DTFS ?
$$
x[n] = \sum_{k=0}^{N-1} a_k e^{jl\omega_o n}
$$

This is, essentially, what the sampling theorem is about.  Uniform sampling in one domain (e.g. the time domain) causes periodic extension in the reciprocal domain (e.g. frequency domain).
The reason why is that the sampling function is a periodic function which means it can be represented as a Fourier series
$$\begin{align}
\operatorname{III}_T(t) \ &\triangleq\  \sum_{k=-\infty}^{\infty} \delta(t - k T) \\
 &= \sum_{n=-\infty}^{+\infty} c_n e^{j 2 \pi n \frac{t}{T}}  \\
\end{align}$$
where the Fourier coefficients are
$$ \begin{align}
c_n\, & = \frac{1}{T} \int_{t_0}^{t_0 + T} \operatorname{III}_T(t) e^{-j 2 \pi n \frac{t}{T}}\, \mathrm{d}t \quad ( -\infty < t_0 < +\infty ) \\[4pt]
& = \frac{1}{T} \int_{-\frac{T}{2}}^{\frac{T}{2}} \operatorname{III}_T(t) e^{-j 2 \pi n \frac{t}{T}}\, \mathrm{d}t \\[4pt]
& = \frac{1}{T} \int_{-\frac{T}{2}}^{\frac{T}{2}} \delta(t) e^{-j 2 \pi n \frac{t}{T}}\, \mathrm{d}t \\
& = \frac{1}{T} e^{-j 2 \pi n \frac{0}{T}} \\[4pt]
& = \frac{1}{T} \ .
\end{align} $$
Using the definition of Fourier Transform most common with electrical engineering
$$ X(f) \triangleq \int\limits_{-\infty}^{+\infty} x(t) e^{-j 2 \pi f t} \ \mathrm{d}t $$
$$ x(t) = \int\limits_{-\infty}^{+\infty} X(f) e^{+j 2 \pi f t} \ \mathrm{d}f $$
when sampling $x(t)$
$$\begin{align}
x_\text{s}(t) &\triangleq x(t) \cdot T \cdot \operatorname{III}_T(t) \\
 &= x(t) \cdot T \sum_{k=-\infty}^{\infty} \delta(t - k T) \\
 &= T \sum_{k=-\infty}^{\infty} x(t) \ \delta(t - k T) \\
 &= T \sum_{k=-\infty}^{\infty} x(kT) \ \delta(t - k T) \\
\end{align}$$
which shows how $x(t)$ is converted to samples $x(kT)$ and this is also true:
$$\begin{align}
x_\text{s}(t) &\triangleq x(t) \cdot T \cdot \operatorname{III}_T(t) \\
 &= x(t) \cdot T \cdot \sum_{n=-\infty}^{+\infty} \frac1T e^{j 2 \pi n \frac{t}{T}} \\
 &= \sum_{n=-\infty}^{+\infty} x(t) \ e^{j 2 \pi n \frac{t}{T}} \\
\end{align}$$
the resulting Fourier transform is
$$ X_\text{s}(f) = \sum_{n=-\infty}^{+\infty} X\left(f - n\tfrac{1}{T} \right)
 $$
which is periodic in the frequency domain with period $\frac{1}{T}$.
Because of the Duality of the Fourier Transform, if can be also shown that sampling in the frequency domain causes periodicity in the time domain, which is essentially all that Fourier series is about.


Transposed FIR filter structure to compensate for extra delays on addition The problem I'm having concerns an add operation in verilog which is part of a transposed FIR filter.

add(a,delay_line[i],delay_line[i+1]);

a is an input from multiplying a cofficient with the current input sample.
delay_line[i] is an input from the delay line at index i
delay_line[i+1] is an output to the delay line at index i+1

The problem I have is that the add operation takes $3$ clock cycles. My new a input arrives every $1$ clock cycle, so as you can see ideally
I would need the add operation to take one clock cycle, then the
inputs to the add are at the same rate, and can update the delay_line
correctly. I'm sure doing an add in one clock cycle is fine in practice,
but in the case where it is not, then how can this be resolved ?.

My first approach was to draw out the filter structure and add more
delay taps on the ouputs of adds. The normal FIR transpose filter has
single $z^{-1}$ delays after each addition, so I placed some $z^{-3}$ delays to compensate for the add delays. doing that wrecks the frequency response so its not that simple :)

I initially thought this was a Verilog / FPGA problem, but I think this
is more of a DSP algorithm problem now so I posted it to this forum. If someone can point me to some literature on this that would be
appreciated.

Here's my attempt to model this in MATLAB (sorry if you don't have this tool)

%this is an fir transpose filter
fre = ( [0 1500 1700 4000]/4000) ;
msk = [1 1 0 0]; %mask
b = firpm(3,fre,msk); %simple 4 coefficient low pass

%convert filter to transposed form
Hd = dfilt.dffirt(b);
b = cell2mat(Hd.coefficients);

bd2 = [b(1) 0 0 b(2) 0 0 b(3) 0 0 b(4)]; %never this simple!

fvtool(b,1,bd2,1);%compare frequency response

Ah, its easier than I thought, just mirror the introduced delays. So this means introduce two additional $z^{-2}$ in the last two stages of the transpose FIR filter structure. After some basic math to get $H(z)$ the coefficients are as below in bd2, which gives the same magnitude frequency response as the original.

bd2 = [0 0 0 0 b(1) b(2) b(3) b(4)];


If there is a function $H(e^{j\omega})$. Then will its conjugate be $H(e^{-j\omega})$?Why? If there is a function $H(e^{j\omega})$. Then will its conjugate be $H(e^{-j\omega})$?Why?
Let there be a function
$$H(e^{j\omega}) = \frac{1+ae^{-j\omega}}{1+\alpha e^{-j\omega}+\beta e^{-j2\omega}}$$
Then will its conjugate be:
$$H^*(e^{j\omega}) = \frac{1+ae^{+j\omega}}{1+\alpha e^{+j\omega}+\beta e^{+j2\omega}}$$
Where $j$ is an imaginary number equal to $\sqrt {-1}$.
How is it valid for the general case?

What you say only holds true for real signals. Let's prove it by the definition of Fourier transform:
$$
H(e^{j\omega}) = \mathcal{F}[h(t)] = \int_{-\infty}^{\infty} h(t) e^{-j\omega t} dt
$$
and $H^*(e^{j\omega})$ should be
$$
H^*(e^{j\omega}) = \int_{-\infty}^{\infty} h^*(t) e^{j\omega t} dt
$$
thus
$$
H^*(e^{-j\omega}) = \int_{-\infty}^{\infty} h^*(t) e^{-j\omega t} dt = \mathcal{F}[h^*(t)]
$$
If $h(n)$ is real, $\mathcal{F}[h^*(t)] = \mathcal{F}[h(t)] $, we can derive that
$$
H^*(e^{-j\omega}) = H(e^{j\omega})
$$
Take conjungate on both sides of the equation above, we have
$$
H(e^{-j\omega}) = H^*(e^{j\omega})
$$


How FFTW factorizes array length? To compute fft with FFTW algorithms I think that it needs to factorizes the array size in prime factors, but it is an hard problem. which Algorithm uses FFTW to do it?

FFTW factorizes by testing against a table of primes.  This is fast due to the fact that the size of the table is limited because the size of computer memory in which the FFT data is to be computed is itself limited to a very tiny number in mathematical terms (less than a few petabits). You only need up to the square root of that number.


The condition in order to not lose information after upsampling and downsampling If you have a signal x[n] and it is interpolated by N and then decimated by M, how do you ensure minimum distortion in x[n] from the resampling process?

The answer to this is to design a filter that perfectly passes the original spectrum and completely rejects the images that the decimation and interpolation process creates. Such a filter is not achievable, so the outcome becomes a design constraint on how much distortion we can tolerate (as in all filter design). Helpful for this is to understand the mechanism that creates the images in this case so I am including some graphics below that I have that further explain this process.

Interpolation by inserting zeros

We upsample with minimum distortion by inserting zeros in between each sample in the time domain, and then follow this with an interpolation filter. The interpolation filter passes the original spectrum and rejects the images the zero-insert creates.

The graphic below shows the image replication of the spectrum by upsampling by 4 (insert 3 zeros between each sample). Note how we started off with just our desired spectrum in the "first Nyquist zone" which is the freqency from DC to half the sampling rate ($F_s/2$) for a real signal. (The same applies to the spectrum for a complex signal, just in that case the spectrum is unique from $F_s/2$ to $+F_s/2$ so we would include that or equivalently the signal from DC to $F_s$.) After the zero insert we end up with 4 additional spectral images in the new first Nyquist zone at the higher sampling rate, for the case of a 4x sampling rate increase. (So in your case you would have N such images).

The ideal interpolation filter would "grow" the zeros to the actual value needed for perfect interpolation. Such a filter could be a multiband filter (Matlab firls and firpm do these easily) as depicted in the graphic below, but a low pass filter would also be fine to use. The actual distortion depends on how much of the image frequencies we allow through and how much pass band distortion (either loss or ripple) we allow, this all becomes a design constraint on the complexity of the filter, tolerable delay etc...

Decimation by dropping samples

For decimation we reduce the sampling rate by simply removing samples (So if we decimate by M, we only use every Mth sample). The spectral challenge of images in the decimation process is very similar. While in interpolation we were concerned about clean-up in that we filter after we do the zero-insert; for decimation we are concerned with filtering prior to the dropping of samples due to aliasing; signal energy that will fold into our spectrum of interest if we don't filter prior to removing samples.  The graphic below shows this for an example of decimating by 4. This is identical to the aliasing issues in A/D conversion (in fact A/D conversion is effectively decimating from an infinite sampling rate! So here we are doing a digital to digital conversion):

Like the interpolation filter, the ideal decimation filter would pass our spectrum of interest with no distortion while completely eliminating any signal that may exist in the higher spectral locations that can fold in. Note if you know by design that these spectral locations have already been filtered (including amplified noise floor), then filtering may not be necessary.

Like the interpolation filter design, a multiband filter approach would provide more rejection for the same number of taps (more efficient):

And finally a rational rate conversion (both decimate and interpolate) suggests that you can do the interpolation first so that the filtering can be shared between the two operations. This is a consideration dependent on the actual filtering design requirements for both operations, but feasible that one low pass filter design could cover both requirements. In some cases with high numbers for M and N it may also makes more sense to break up the interpolation and decimation operations further to simplify the filtering requirements.

In all cases a significant limit to what can be achieved is how much of the digital spectrum available we are using. The more excess spectrum available, the simpler the filter design becomes since the design is driven by how tight the transition bands need to be.


How to evaluate quality of signal? I'm trying to evaluate received signals from a transceiver. I want to make a threshold of received signals so that we can say some signals are meaningless.

I have learnt something about Signal-to-Noise-Ratio (SNR) and Received-Signal-Strength (RSS). But I found that symbols can still be decoded even in low SNR ($\le 0$ dB).

Here are my questions:

If I choose SNR as my metrics, will my threshold be relevant to modulation or coding schemes?
Will RSS be suitable for theoretic analysis as metrics?

Indeed, in practical systems there is usually a cut-off SNR above which communication is possible with a given amount of reliability. This is usually determined according to a design goal, expressed in terms of a bit error rate (BER) or a packet error rate (PER). Typically, the BER should not exceed $10^{-4}$ or $10^{-5}$.
According to your signal parameters (modulation and coding), sometimes it is possible to derive the BER as a function of the SNR. This depends on the complexity of your transmitter-receiver architecture and the channel model you use.
So:

Yes.

Now, regarding the RSS, I would say that it is not a standard metric. Some equipments give RSS metrics, according to an unknown criteria. Sometimes, it is just a table where a given range of received power is indexed. Hence:

It depends. In general, we don't use the RSS to directly determine design parameters of a demodulator (in this case, a threshold). But if you know exactly what it means, you may derive a relation SNR-RSS. But since RSS measures the total received signal power--including noise--I suppose you will have to estimate the noise power also.


Sampling Frequency and Spectral Regrowth Sampling a cosine wave of 10 Hz at Fs = 64 and number of samples Ns = 256
I setup my time vector for the cosine wave as n=(0:Ns-1)*(1/Fs)

If I change the sampling frequency from 64 to 64.0005 I get imaginary spectrum coming up from the quantization noise and actually having significant amplitudes. I also get this at 65, 66, 67, etc but not at powers of 2, ie Fs=128,256...

I thought it was because the signal was truncated, but I substantially increased the number of samples (Ns) and this had little to no effect in removing the imaginary components of the cosine.

Can someone please advise what the link here is?

Yes that's called as the aperture effect and indicates that if the FFT window length (the aperture) coincides with an integer number of signal periods, then the resulting FFT spectrum will have no leakage artefacts (mostly having zeros, up to quantization noise, bins) and significantly nonzero bins only at the signal frequencies.


Choosing the best measure for machine learning on a Spectrogram I am working on an application to classify Talas(a piece of music with onsets). For that, I have the spectrogram of a Tala, which I want to compare with the test Tala, or compare 2 Talas.

For this I intend to find the onsets in the spectrogram and then form a classifier based on the frequencies that I get at the onset points in the Tala.

Now as we already know that a spectrogram contains multiple frequencies for a time t , I wish to quantify this as a single value to be used in machine learning, so that I have n values for the n onsets(say), and not multiple frequency values for the same time t.

So here is my question - Given all the frequencies at time t, how can I get a value that represents all these frequencies(like mean, median, etc). What should be the measure I should go ahead with?

First thing: a spectrogram does not give the "frequencies" at time $t$. It gives amplitudes, measured for a set of linearly spaced frequency bands.

Anyway, there is no "right" answer to your question, because it depends on which musicological aspect you want to focus on, and you need to enlighten us on that. If you want to focus on melody, extracting the fundamental frequency is the most relevant feature, but I suspect this is not the case here. If you want to discriminate different timbres, spectral centroid might be a better way of "summarizing" a slice of spectrogram into a single value - better, but still quite bad.

Indeed I can see several flaws in your methodology:

Why extract information at time t only (assumedly at a small window centered at t) rather than on the entire segment of sound between onsets?
Why do you have this constraint that for any onset you have to extract a single real value? Some properties of sound, like timbre, are multidimensional, so you won't be able to entirely capture them in one single value. If we take the example of discriminating tabla sounds (which might or might not relevant to your problem, since you mentioned tala...), we need several features to discriminate them (bayan/dayan discrimination relies on pitch information, while resonant vs non-resonant stroke discrimination relies on timbral features).
You are completely discarding the rhythmic (time) information. Is that really wanted for your problem?

Tell us more about your problem, there's a good chance it has already been solved, and that you are asking the wrong questions to solve it...


Estimate digital signal I have a device that connect with USB to my computer. This device have digital signal that I can sample it every few seconds. The signal that the device perform, is a digital signal with $X$ seconds of $1$ and $Y$ seconds of $0$ ($X$ and $Y$ are integers and consts).

What I want to do is to estimate if in $Z$ seconds from now, the signal will be $1$. The problem is that I don't know what are $X$ and $Y$ (I don't know how much time $1$ and $0$ are), so I don't know how can I estimate it, and I can sample it every 10 seconds ($X$ and $Y$ can be less then 5).
But if I sample $0$ at some moment, the device tells me in how much seconds $1$ will appear again.

I have basic knowledge in signal analysis and I don't know how to do that.

Is it possible at all? If it is, how can I do it?

If X plus Y is longer than twice your sample rate, 20 seconds, but not an exact rational multiple of your sample rate, then with a large enough number of samples you might be able to make a frequency and duty cycle estimate by calculating the statistics and spectrum of the sample vector.


Difference between step,ramp and Impulse response As, we can determine the response of the system by its 'impulse response' then why we use 'step response' and 'ramp response'. Is there any difference between all these response?

Like the impulse response, the step and ramp responses describe an LTI-system. There are situations when the step response or the ramp response describe what we want to know about the system more directly than the impulse response. When the input to a system is a piece-wise constant function of time, the transitions are well described by the step response. It gives quite directly information about overshoot and settling time.


Meaning of Q and I in all-software QPSK communication systems? Q and I seem to come up a lot in discussions of digital communications, such as the comprehensive answer provided to this question:

FFT window phase syncronising

It seems like the Q and I channels ultimately are combined into a sin wave of a given frequency, phase and amplitude.  If the communications system is all software (for example in a software-defined radio) is it necessary to worry about "Q" and "I"?  Is this just an artifact of when receivers and transmitters were electronics based?

To make the question more concrete, the link provided talks about putting the even bits onto Q and the odd bits onto I.  If this is the same thing as creating a sin wave of phase 0, 90, 180, or 270, then why bother with this seemingly unnecessary notation?  Perhaps the pulse shaping (like an RRC filter) needs to be independently applied to the Q and I channels?

A point in 2D space can be represented by either X and Y coordinates, or a distance and angle (relative to some reference point and axis).  However, the point measurement errors or noise is some common real systems is more likely to be uniform in the the units of an X,Y coordinate system, which may make error analysis and numerical error control in this space more tractable, thus lowering the computational requirement for most types of processing, even in purely digital systems.


Maximum Likelihood for Colored Noise I have the following question about the maximum likelihood (ML) in presence of inter-symbol interference and colored noise.

Assume the communication system is as follows. Information source, modulator, transmit pulse filtering, channel, AWGN, matched receiver pulse filtering, sampler, and then ML estimator.

The received signal after the sampler is given as

$$y = h\star h\star c\:a + n\star h$$

where $y$ is the received samples (contains transmit symbols + ISI + noise), $a$ is the transmit symbols, $h$ transmit/receive pulse shaping filter, $c$ is the channel impulse response, $n$ is the AWGN, and $\star$ denotes convolution. This can be further written as

$$y = A\:a + \eta$$

where $A$ is a convolution matrix whose elements contains $h\star h\star c$ and $\eta$ is a colored noise.

I wonder what will be the ML in this case?
And how to derive/calculate it?

Thanks.

Let's have a look on the following model:

$$ y \left[ n \right] = \left( h \ast x \right) \left[ n \right] + \left( g \ast w \right) \left[ n \right] $$

Where $ x \left[ n \right] $ is the signal of interest and $ w \left[ n \right] $ is the AWGN with unit Variance.

In Matrix form it is written by:

$$ \boldsymbol{y} = H \boldsymbol{x} + G \boldsymbol{w} $$

Where $ H $ and $ G $ are the convolution matrices of the model.

Let's define $ v = G \boldsymbol{w} $ then $ v \sim \mathcal{N} \left( \boldsymbol{0}, G {G}^{T} \right) $ which implies (If the input to Linear Operator is Gaussian Variable then the output is also Gaussian Variable):

$$ \boldsymbol{y} \sim \mathcal{N} \left( H \boldsymbol{x}, G {G}^{T} \right) $$

Then the Maximum Likelihood is given by:

$$ \arg \max_{ \boldsymbol{x} } \det \left( 2 \pi G {G}^{T} \right)^{-\frac{1}{2}} {e}^{ -\frac{1}{2} {\left( H \boldsymbol{x} - \boldsymbol{y} \right)}^{T} {\left( G {G}^{T} \right)}^{-1} \left( H \boldsymbol{x} - \boldsymbol{y} \right) } $$

It is easy to see you only need to deal with the term in the power of the Exponent, so it is equivalent of:

$$ \arg \min_{ \boldsymbol{x} } {\left( H \boldsymbol{x} - \boldsymbol{y} \right)}^{T} {\left( G {G}^{T} \right)}^{-1} \left( H \boldsymbol{x} - \boldsymbol{y} \right) $$

Which is nothing more than Weighted Least Squares problem which matches the result of Linear Regression for Colored Noise (As expected for Gaussian Noise).

The optimal solution is given by:

$$ \hat{x} = {\left( {A}^{T} {\left( G {G}^{T} \right)}^{-1} A \right)}^{-1} {A}^{T} H b $$


noise cancellation of non-repetitive noise There are a new noise cancellation headphones that almost fully cancel human noise.

What kinds of algorithms could achieve that ?

EDIT:by human noise , i mean talking. In general, previous audio cancellation headphones didn't work that well for attenuating talking.

Since I cannot comment, here comes a short answer.

As far as I know, there are microphones outside the earphone. The microphones record the outside noise and generate sounds with 180 phase shift to cancel the noises. The noises are canceled actively.


Fourier Transform of Kernel Density Estimation - Convolution Theorem? I am reading this paper about density estimation (Appendix A), where the authors apply a Fourier transform to the estimated probability density (the $X_j$ are a sample of $N$ data points drawn independently from an unknown probability density function $f(x)$):

$$
\hat{f}(x) = \frac1N \sum^N_{j=1} K(x-X_j)
$$

The authors now note that they apply the convolution theorem and the Fourier transform (where the Fourier transform of the function $g(x)$ is defined by $\Phi_g(t) = \int \exp\lbrace\mathrm{i}xt\rbrace g(x) \mathrm{d}x$), and obtain:

$$
\Phi_{\hat{f}}(t) = \kappa(t) \frac1N \sum^N_{t=1} \exp\lbrace\mathrm{i}tX_j\rbrace
$$

where $\kappa(t)$ is the Fourier transform of the Kernel $K$.

Now, I am a bit stuck on how they obtained this form and how exactly they used the convolution theorem here. In particular, what happened to the argument $X_j$ in the kernel?

EDIT: Ok, I think I got it. Is this the right approach?

We apply the Fourier transform and obtain:

$$
\Phi_{\hat{f}}(t) = \frac1N \sum^N_{j=1} \int K(x-X_j) \exp\lbrace\mathrm{i}xt\rbrace\mathrm{d}x
$$

We can easily rewrite this equation as an integral over a Dirac delta distribution:

$$
\Phi_{\hat{f}}(t) = \frac1N \sum^N_{j=1} \iint K(y) \delta(y-(x-X_j)) \exp\lbrace\mathrm{i}xt\rbrace\mathrm{d}x\mathrm{d}y
$$

Integrating over $x$ and using the symmetry of $\delta$ we obtain:

$$
\Phi_{\hat{f}}(t) = \frac1N \sum^N_{j=1} \int K(y) \exp\lbrace\mathrm{i}(y+X_j)t\rbrace\mathrm{d}y
$$

We see that we can move factor involving $X_j$ outside the integral, and are left with the result:

$$
\Phi_{\hat{f}}(t) = \frac1N \sum^N_{j=1} \exp\lbrace\mathrm{i}X_jt\rbrace \int K(y) \exp\lbrace\mathrm{i}yt\rbrace\mathrm{d}y
$$

where the first factor is $\Delta(t) = \frac1N \sum^N_{j=1} \exp\lbrace\mathrm{i}X_jt\rbrace$, and the second factor is $\kappa(t) = \int K(y) \exp\lbrace\mathrm{i}yt\rbrace\mathrm{d}y$, exactly as in the paper.

I am still not sure in how far they applied the convolution theorem here?

It's a bit contrived, but observe that by the "sifting property" of Dirac-delta function:

$$
K(x-X_j) = K(x) * \delta(x-X_j)
$$

where $*$ denotes convolution and $\delta$ is the Dirac-delta function.

Now you can apply Fourier convolution theorem and write:

$$
\mathcal{F}\{K(x-X_j)\}(t) = \kappa(t) e^{j t X_j}
$$

where $\mathcal{F}\{g(x)\}(t)$ denotes the Fourier transform of a function $g(x)$.

Finally using the linearity property of the Fourier transform:

$$
\mathcal{F}\{\hat{f}(x)\}(t) = \frac{1}{N} \sum_{j=1}^N \mathcal{F}\{K(x-X_j)\}(t) = \kappa(t) \sum_{j=1}^N e^{j t X_j}.
$$


How is the simplified version of the Bromwich inverse Laplace transform integral derived? I do not understand how the last equality is derived from the previous.
Apparently the first term in the integral (involving $\mathrm{cos}$) is equivalent to the second (involving $\mathrm{sin}$)!! How so??

I DO understand how the integral range is halved (since $F(s)^*=F(s^*)$; where $F(s)$ is the Laplace transform of $f(t)$.
Any help would be appreciated since this form is used often in numerical inverse Laplace transform algorithms.

[Note: $\hat{f}(s)$ below represents the Laplace transform of $f(t)$]

This quote is from the web source Abate and Whitt, 1995.

I agree that the derivation is unclear, yet the final result is correct (for $t>0$, see below). There are two conditions that are necessary for the final result to be true:

$f(t)$ is real-valued
$f(t)$ is causal

The step from line 2 to line 3 in the derivation assumes that $f(t)$ is real-valued, i.e., only the real part of the integrand is considered. The last step leading to the final result assumes causality of $f(t)$, i.e., $f(t)=0$ for $t<0$.

From the third line in the derivation we have for real-valued $f(t)$

$$f(t)=\frac{e^{at}}{2\pi}\int_{-\infty}^{\infty}\left[\text{Re}\left(\hat{f}(a+iu)\right)\cos(ut)-\text{Im}\left(\hat{f}(a+iu)\right)\sin(ut)\right]du\tag{1}$$

and, consequently,

$$f(-t)e^{2at}=\frac{e^{at}}{2\pi}\int_{-\infty}^{\infty}\left[\text{Re}\left(\hat{f}(a+iu)\right)\cos(ut)+\text{Im}\left(\hat{f}(a+iu)\right)\sin(ut)\right]du\tag{2}$$

Since for causal $f(t)$ we have $f(-t)=0$ for $t>0$, we can write

$$f(t)=f(t)+f(-t)e^{2at},\qquad t>0\tag{3}$$

Consequently, adding $(1)$ and $(2)$ gives

$$f(t)=\frac{e^{at}}{\pi}\int_{-\infty}^{\infty}\text{Re}\left(\hat{f}(a+iu)\right)\cos(ut)du,\qquad t>0\tag{4}$$

And since for real-valued $f(t)$ the integrand in $(4)$ is even, we finally obtain

$$f(t)=\frac{2e^{at}}{\pi}\int_{0}^{\infty}\text{Re}\left(\hat{f}(a+iu)\right)\cos(ut)du,\qquad t>0\tag{5}$$

q.e.d.

Note that this is result is only valid for $t>0$, which is not stated in the paper you quoted. Of course, for $t<0$ we have $f(t)=0$.

Also note that instead of $(3)$ we could have written

$$f(t)=f(t)-f(-t)e^{2at},\qquad t>0\tag{6}$$

from which we can conclude that

$$f(t)=-\frac{e^{at}}{\pi}\int_{-\infty}^{\infty}\text{Im}\left(\hat{f}(a+iu)\right)\sin(ut)du,\qquad t>0\tag{7}$$

and, taking into account that for real-valued $f(t)$ the integrand is even, we get

$$f(t)=-\frac{2e^{at}}{\pi}\int_{0}^{\infty}\text{Im}\left(\hat{f}(a+iu)\right)\sin(ut)du,\qquad t>0\tag{8}$$

Comparing $(5)$ with $(8)$ we see the equivalence

$$\int_{0}^{\infty}\text{Re}\left(\hat{f}(a+iu)\right)\cos(ut)du=-\int_{0}^{\infty}\text{Im}\left(\hat{f}(a+iu)\right)\sin(ut)du\qquad t>0\tag{9}$$


Help with obtaining the power spectral density of a simple continuous cosine (using both forms of the definition for PSD) I am trying to go through a simple example to teach myself about Parseval's theorem and calculating power spectral density (PSD) in practice and would be very grateful if someone could check my reasoning and help my understanding.

Specifically, I want to calculate the average power of a signal in the time domain and show that it is equal to the average power obtained in the frequency domain using the PSD (according to Parseval).

As an example, I am considering a simple cosine (non-causal) signal $x(t) = A\cos(2\pi f_0t)$, which should have infinite energy but finite average power (known as a "power signal", as opposed to "energy signal") given by:
$$P_{\textrm{av}} = \lim_{T\to\infty}\frac{1}{T} \int^{+T/2}_{-T/2} |x(t)|^2\mathrm dt$$

Since this signal is periodic, I should be able to calculate the average power by considering a single period only, where $T= 1/f_0$,
$$P_{\textrm{av}} = \frac{1}{T} \int^{+T/2}_{-T/2} |A\cos(2\pi f_0t)|^2\mathrm dt = f_0 A^2 \int^{+T/2}_{-T/2} \frac{1}{2}\Big[1+\cos(4\pi f_0 t) \Big]\mathrm dt = \frac{A^2}{2}$$

I would now like to arrive at this result by integrating the power spectral density over all frequencies (as should work by Parseval), to convince myself of what I'm doing. So first, I need to obtain the power spectral density. I have seen one definition of the PSD given as the Fourier transform of the autocorrelation function, $R(\tau)$, so I first calculate this:

\begin{align}
  R(\tau) &= \int^{+\infty}_{-\infty} x(t+\tau)\;x^*(t)\;\mathrm dt  \\
          &= A^2 \int_{-\infty}^{+\infty} \cos(2\pi f_0(t+\tau))\cdot \cos(2\pi f_0)\; \mathrm dt\\
          &= \frac{A^2}{2} \cos(2\pi f_0\tau)
\end{align}

where I have used trigonometric identity to evaluate the integrals. Now, calculating the Fourier transform of this to get the PSD:

\begin{align}
  \textrm{PSD}(f) &= \mathcal{F}\{R(\tau)\}   \\
                  &= \int_{-\infty}^{+\infty} R(\tau) e^{-2\pi i f \tau}\; \mathrm d\tau\\
                  &= \int_{-\infty}^{+\infty} \frac{A^2}{2} \cos(2\pi f_0\tau)  e^{-2\pi i f \tau}\; \mathrm d\tau\\
                  &= \frac{A^2}{4}\Big[ \delta(f-f_0) + \delta(f+f_0) \Big]
\end{align}

Is this correct for the power spectral density of a cosine wave, i.e. in units of [signal$^2$ per Hz]? It does indeed look like if I were to integrate this PSD over frequency I would get the correct average power $P_\textrm{av} = A^2/2$.

I have seen an alternative (or just different form?) of the definition of PSD in this question:

$$S_{xx}(\omega)=\lim\limits_{T\to \infty}\mathbf{E} \left[ | \hat{x}_T(\omega) |^2 \right]$$

How would I apply this definition to my cosine signal to arrive at the same PSD above, and show that the average power is recovered? Which method is the approach I should take? Is it true that the autocorrelation method is used more for stochastic signals when the FT does not exist, and for deterministic signals (such as in my case) we can directly use the FT?

There are several misconceptions in the question that have not been addressed in the existing answers. First of all, the signal $x(t)=A\cos(2\pi f_0t)$ is a deterministic power signal (unless $A$ or $f_0$ are modeled as random variables). For this reason several definitions in the question are inappropriate. First, the auto-correlation of a power signal is given by

$$R_x(\tau)=\lim_{T\to\infty}\frac{1}{2T}\int_{-T}^Tx^*(t)x(t+\tau)dt\tag{1}$$

The integral given in the question (with infinite limits and without division by $T$) does not exist for the given $x(t)$. With definition $(1)$, the auto-correlation of $x(t)$ is indeed obtained as

$$R_x(\tau)=\frac{A^2}{2}\cos(2\pi f_0\tau)\tag{2}$$

The Fourier transform of $(2)$ results in the power spectrum of $x(t)$.

The power spectrum can also be computed directly from $x(t)$, but the formula given in the question only applies to random signals, but not to deterministic signals. For deterministic signals, the appropriate definition is

$$S_x(f)=\lim_{T\to\infty}\frac{1}{2T}\left|\int_{-T}^{T}x(t)e^{-j2\pi ft}dt\right|^2\tag{3}$$

The computation of $(3)$ for the given signal is discussed in this question.


Automatic Image Enhancement of Images of Scanned Documents (Auto Whitening) Dropbox have make a blog post about there automatic enhancement method for scanned document image - Fast Document Rectification and Enhancement.

I followed the post and they mention a formula to make the enhancement:

J(x,y) = gain(x,y)*I(x,y) + offset(x,y)

Where $ J $ is the output image and I is the Input image, the problem is i don't understand how they calculate gain and offset for each pixel

Thank you for your reply

EDIT: add simple image

Their intention is to solve:

$$ \arg \min_{ J \left( x, y \right) } \sum_{x, y} \left\| J \left( x, y \right) - 255 \right\|_{2}^{2} + \lambda \sum_{x, y} \left\| \nabla J \left( x, y \right) - \nabla I \left( x, y \right) \right\|_{2}^{2} $$

Where $ I \left( x, y \right) $ is the input image and $ J \left( x, y \right) $ is the output image.

They state that usually the solution is constant locally and can be approximated by the map of Gain and Offset.
They do not say the method they use to calculate it but you can use few ideas:

Solve the above optimization problem for a downsampled version of $ I \left( x, y \right) $. From the downsampled solution of $ J \left( x, y \right) $ infer the downsampled map of the Gain and Offset. Then upsample the maps and apply them on the full size $ I \left( x, y \right) $.
Just like the Guided Filter approximate a local solution of the Weighted Least Squares filter, reformulate it for solution of the above problem (You'll see they are very similar).

Solving the Optimization Problem

As written in Dropbox's Blog Post - Fast Document Rectification and Enhancement the above is a Poisson Equation which has many dedicated solvers.
To give you a glimpse of straight forward way solve it we can vectorize the image into a vector. Let's say the input image is $ y $ and the output image is $ x $ (Vectors) then the problem can be rewritten:

$$ \hat{x} = \arg \min_{x} \frac{1}{2} {\left\| x - 255 \cdot \boldsymbol{1} \right\|}_{2}^{2} + \frac{\lambda}{2} {\left\| {D}_{h} \left( x - y \right) \right\|}_{2}^{2} + \frac{\lambda}{2} {\left\| {D}_{v} \left( x - y \right) \right\|}_{2}^{2} $$

Where $ {D}_{h} $ is the horizontal Derivative Operator, $ {D}_{v} $ is the vertical Derivative Operator and $ \boldsymbol{1} $ is vector of ones.
Then the solution is given by:

$$ \hat{x} = { \left( I + \lambda {D}_{h}^{T} {D}_{h} + \lambda {D}_{v}^{T} {D}_{v} \right) }^{-1} \left( \lambda {D}_{h}^{T} {D}_{h} y + \lambda {D}_{v}^{T} {D}_{v} y + 255 \cdot \boldsymbol{1} \right) $$

Pay attention that if the image $ I \left( x, y \right) \in \mathbb{R}^{m \times n} $ then $ D \in \mathbb{R}^{m n \times m n} $.
Hence from memory and computation this is hard to do for large images unless the solver takes advantage of the matrix $ D $ structure (Sparse) and that's what Poisson Solvers are all about (Be more accurate $ {D}^{T} D $ has a structure to take advantage of without ever formally creating it).
For very small images you can play with this.
For large you can do follwoing:

Use Direct Solvers
For example just using the \ operator in MATLAB.
Use Iterative Solvers
For example the Conjugate Gradient. Nice property of the Conjugate Gradient is it doesn't require the matrices themselves but the products of the matrices with the image which can be calculated by convolution.

So basically the above is about solving Linear Equation with special structure. You may read about solving those kind of equations at the course The Poisson Equation in Image & Shape Processing.

MATLAB Code

I added a MATLAB Code to solve the problem of whitening the image as described in the Blog Post.

Here is the result for different values of $ \lambda $:

As expected, one could see that the higher the regularization the better the contrast.

The full MATLAB code which implements the solution is at my Signal Processing StackExchange Question 50329 - GitHub Repository.
The code uses direct solver. In practice one should use Preconditioned Conjugate Gradient solver which only requires the result of the operator on the image which can be done using Convolution in much faster manner than Matrix Multiplication (At least for "Denser" operators).

As reference image I used Certificate of Arrival for Berta Werner.


Presenting an all-pass filter of a high order as a cascade of all-pass filters of the 1-st order Are all-pass filter of the n-th order just a cascade of N all-pass filters of the 1-st order?

There are some nice statements regarding zeros, poles, and gain for the 1-st order. Can those statements be applied to the n-th order all-pass filter too? For instance, the Pole is a reverse complex conjugate of the Zero.

Obviously, cascade of N all-pass filters of the 1-st order will produce the all-pass filter of n-th order, so the 'necessary' condition is true. But is the 'sufficient' condition true as well?

There are some edge cases of course, for example, if we consider two kind of all-pass filters with abs(H(exp(iw)))=1/2 and abs(H(exp(iw)))=2 then their multiplication (i.e. cascade) will give a true all-pass filter. But to me it's still an edge case. I'm more interested in non-constant Amplitude Frequency Responses.

This depends a bit how rigorous you define "allpass" filter.

You can show that any pole can be turned into an allpass filter if, and ONLY if, it you pair it with a zero at the inverse location. A zero at the inverse location is the only way to achieve $|H(\omega)|^2 = 1$ for all $\omega$.

Poles can be complex or real. Complex poles result in second order allpass filters, real poles in first order allpass filters.

Pure delays are simply a special case of a single real pole  at z=0 and the matching zero at z=infinity. The transfer function of a one sample delay is almost identical to that of an allpass filter with a pole at .00001 or -.00001. Multi sample delays are simply cascades of single sample delays, so an N-sample delay has N poles at z=0 and N zeros at z=infinity

All of these taken together implies that, indeed, you can represent every allpass filter as a cascade of first and second order allpass filters.

There is a downside to that as well: For first order section the phase is zero at DC and $\pi$ at Nyquist. For a second order it starts at 0 and ends at $2\cdot \pi$ and the phase decreases monotonically between those values. These properties are maintained in a cascade so we can conclude that any N-th order allpass has a phase of 0 at DC and a phase of $N \cdot \pi$ at Nyquist with a monotonically decreasing phase in between. That's too restrictive for many applications. (I'm ignoring a potential multiply with -1 flip here, which would just result on an phase offset of $\pi$).

There is a different class of filters that are "almost" allpass filters, i.e. you can make the magnitude as close to unity as needed over the frequency range of interest. These can not be represented as cascades of 1st and 2nd order allpass filters.

An interesting example of this is a Hilbert transformer: it has unity amplitude gain but a phase shift that's only 90 degrees at Nqyuist. You could interpret this as having an infinite amount of poles and zeroes but that's not very useful. You can't implement an ideal Hilbert transformer, but you most certainly can do something that's "good enough" for your specific application.


Significance of steps involved in flipping process of an image? The above image on the right was obtained by following steps

Multiplying the image on the left by $(-1) ^ {(x+y)}$
Computing the DFT
Taking the complex conjugate of the transform
Computing the inverse DFT
Multiplying the image on the left by $(-1) ^{ (x+y)}$

What is the significance of each of the above 5 steps? Here by significance, I mean: what is the benefit of that step? And what loss can occur if we skip that step?

Please note that Steps 1 and 5  involve a multiplication by $-1$ and that $-1$ is raised to the power $(x+y)$.

Looking at the effect you have achieved; namely a flip along the horizontal and vertical axes, all you need is to apply steps 2,3, and 4.

The following MATLAB / OCTAVE code would also generate the same flipped image in 3 steps as would given by your 5 step procedure:

% I is the image data
If = fft2(I);             % S1 - Take DFT2 of the real image
Gf = conj( If ) ;         % S2 - Take conjugate of DFT of image
g  = real( ifft2(Gf));    % S3 - Invert back into spatial domain

figure,subplot(2,1,1)
imshow(I);
subplot(2,1,2)
imshow(g);

The flipping result can be seen from the following symmetry property of DTFT / DFT :

$$f[n,m] \leftrightarrow F(\omega_1,\omega_2) \implies f^*[-n,-m] \leftrightarrow F^*(\omega_1,\omega_2) $$

For real data $f^*[-n,-m] = f[-n,-m]$.


Gammatone filter vs 2D FFT for spectograms I am working on speech segeration and each source I am reading , spectograms are being made by filterbanks which uses gammatone filters . I can't find a proper c++ source code which would make me understand process behind gammafilters, completely.

Now I am curious what is its difference than 2d FFT, I mean why all sources prefer that way(gammatone filters), but not taking FFT for a time period ?

An FFT is a linearly spaced frequency filter bank of fixed length.  To characterize speech more like the human ear-brain combination, a log or MEL frequency spaced filter bank, with differing filter response lengths for different frequency/pitch wavelengths, is a better match.  An array of gamma tone filter kernels is more like the latter than is an FFT.


Why is the dB plot of a signal usually non-positive? Almost every decibel plot of a signal is usually below 0. Why is this the case? I recently just plotted a transfer function with a very high decibel plot above 0.

In power spectra of signals, 0 dB is some agreed power level. It can for example represent the level of a full-scale (FS) sinusoid. In that case it would be rare to see values above 0 dBFS (sinusoid), as it would require peak cancellation as with $1.125\sin(x) + 0.125\sin(3x)$ to stay within a numerical amplitude range of $[-1, 1]$.  When a power spectrum is given in dB, the 0 dB reference should be specified. Otherwise one can't be sure what it means.

You saw the above 0 dB values not with a power spectrum of a signal but with a transfer function, or equivalently with a magnitude frequency response. In a magnitude frequency response, 0 dB means a magnitude or gain of 1, that there is neither attenuation or amplification. If you followed this standard convention and saw values above 0 dB, it means that the corresponding frequencies are amplified.


How can a system be unstable if $L(j\omega)$ is never exactly $-1$? Say we have a linear system with unity feedback, with loop transfer function $L(j\omega)$. The closed-loop transfer function from reference to output is $T(j\omega) = \frac{Y(j\omega)}{R(\omega)}=\frac{L(j\omega)}{1+L(j\omega)}$.

At frequencies for which $L(j\omega)$ approaches $-1$, clearly $|T(j\omega)| \rightarrow \infty$, so the system is unstable - if excited at this frequency, the output is unbounded.

But systems can be unstable even if $L(j\omega)$ never equals $-1$. From the Nyquist plot, we can  show that $T(j\omega)$ can have poles in the right half plane if the contour encircles $-1$, even if $L(j\omega)$ never exactly reaches it. (However, I have very little intuition for why this is true - I just see it as a theorem from complex analysis that happens to be useful here).

Alternatively, from the Bode plot, we say a system is unstable if there are any frequencies for which $|L(j\omega)|$ > 1 and $\angle L(j\omega)< -\pi$ (i.e. phase margin is negative, or gain margin is less than unity). However, I'm not sure why these two conditions result in instability, since these don't result in $|T(j\omega)|$ going to $\infty$.

Two questions:

(1) If $|T(j\omega)|$ never goes to infinity (which is the case when $L(j\omega)$ is never exactly $-1$), how can a system possibly be unstable? Is $|T(j\omega)| \rightarrow \infty$ not the right criterion for deciding whether a system is unstable?

(2) Intuitively, why is a system unstable if there are any frequencies for which $|L(j\omega)|$ > 1 and $\angle L(j\omega)< -\pi$? I understand that you can see it from the Nyquist plot because these two conditions tend to result in encirclements of $-1$ in the $L(s)$ plane, but I'm looking for the intuition.

You cannot make conclusions about the stability of a system by only considering its transfer function evaluated on the imaginary axis $s=j\omega$. Replacing $s$ by $j\omega$ in the transfer function only makes sense for a stable system, otherwise you get a function of $\omega$ that does not describe the system, but another (stable) system.

Let me explain this by an example. Assume a causal system with a transfer function

$$H(s)=\frac{1}{s+a}\tag{1}$$

If $a>0$, the system is stable and by substituting $s=j\omega$, we get its frequency response

$$H(j\omega)=\frac{1}{j\omega+a}\tag{2}$$

The corresponding impulse response is

$$h(t)=e^{-at}u(t)\tag{3}$$

where $u(t)$ is the unit step response. We have

$$H(s)=\int_{-\infty}^{\infty}h(t)e^{-st}dt,\qquad \text{Re}\{s\}>-a\tag{4}$$

and

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

Now imagine that we decrease the value of $a$. The system remains stable as long as $a>0$, but as soon as $a$ becomes $0$, $H(j\omega)\rightarrow\infty$ for $\omega=0$, and the system becomes unstable. If we further decrease $a$, the system remains unstable, but the expression $(2)$ does not go to infinity anywhere. However, the point is that the expression $(2)$ is not the frequency response of the system anymore if $a<0$. The reason is that the integral $(5)$ does not exist anymore. So it is pointless to try to deduce stability or instability from the frequency response of a system, because an unstable system has no frequency response.

What you get if you replace $s$ by $j\omega$ in the transfer function of an unstable system is the frequency response of another system with the same expression for the transfer function, but with a different region of convergence (ROC). The ROC is the region in the complex $s$-plane for which the integral $(4)$ converges. For the causal system in our example, the ROC was given by $\text{Re}\{s\}>-a$ (the ROC is always a right half-plane for causal systems). For an unstable system with $a<0$, the ROC does not include the $j\omega$-axis, and, consequently, $H(j\omega)$ does not exist. In that case, the expression $(2)$ is the frequency response of a different system with transfer function $H(s)$ given by $(1)$ but with a different ROC given by $\text{Re}\{s\}<-a$. That system is stable but anti-causal, and its impulse response is given by

$$h(t)=-e^{-at}u(-t)\tag{6}$$

which is clearly different from $(3)$.

Concerning your second question, note that $e^{\pm j\pi}=-1$, so a phase shift of $\pi$ for certain frequencies implies a sign change, which turns the negative feedback into a positive feedback. This can be compared to the positive feedback in an audio system where a microphone picks up the loudspeaker signal, which is fed back into the amplifier, sent to the loudspeaker, etc.


Discrete State Space Model - Why Are We Calculating $ x \left[ k + 1 \right] $ Instead of $ \dot{\boldsymbol{x}} \left( t \right) $? A continuous state space model is defined as follows.
$$
\dot{\boldsymbol{x}}(t) = A \boldsymbol{x}(t)+ B \boldsymbol{u}(t) \\
\boldsymbol{y}(t)= C \boldsymbol{x}(t)+ D \boldsymbol{u}(t)
$$

If we want to calculate the state $\boldsymbol{x}$ for a specific time $t$, we need to use the state transition matrix $\Phi$. However, when looking at a discrete state space model, things look different.
$$
\boldsymbol{x}[k+1] = {A}_{d} \boldsymbol{x}[k] + {B}_{d} \boldsymbol{u}[k] \\
\boldsymbol{y}[k] = C \boldsymbol{x}[k] + D \boldsymbol{u}[k]
$$

With the discrete state space model, there seem to be no 1st order differential equations incorporated anymore. Instead, we can directly calculate the next sample, without the need for a state transition matrix. How come this is the case?

I will ask you something that will give you intuition.
How would you calculate the Gradient of an image?
Image is a discretization of reality, so how would you estimate the gradient of the "Reality" if you're given only the image?

In the case above we use Finite Differences to approximate the continuous derivative.
So what actually is approximating $ \dot{x} \left( t \right) $ is:

$$ \dot{x} \left( k \Delta t\right) \approx \frac{ x \left[ k + 1 \right] - x \left[ k \right]}{ \Delta t} $$

I do a little abuse of the $ k $ parameter yet I think you can understand what is going on there.


Is uniform noise present in some real processes? Is noise with uniform distribution present in some real-life applications?

It seems that most of the internet sources describe uniform noise as something theoretical.

Typically, when a real value $x$ from a sensor is discretized ($x_Q$) with a uniform quantizer, the error $x-x_Q$ is likely to be uniformly distributed.

Caveat: as commented by @bone, this can be stated in precise terms, see for instance:

A necessary and sufficient condition for quantization errors to be uniform and white,A. Sripad and D. Snyder, 1977

  In this paper, a necessary and sufficient condition is given to model
  the output of a quantizer as an infinite-precision input and an
  additive, uniform, white noise. The statistical properties of the
  quantization error are studied, and a detailed analysis for Gaussian
  distributed inputs is given.

Statistical Theory of Quantization, Bernard Widrow, Istvan Kollar,  Ming-Chang Liu, 1996

  The effect of uniform quantization can often be modeled by an additive
  noise that is uniformly distributed, uncorrelated with the input
  signal, and has a white spectrum. This paper surveys the theory behind
  this model, and discusses the conditions of its validity. The
  application of the model to floating-point quantization is
  demonstrated.


Recover myself from Dilate and Erode for noise reduction i have obtained some skin detection map (i post below two example of detection map). but, i'm using a normal webcam and the map is noisy. Not a lot, just a bit.
I think i can remove noise only with some combination of dilate and erode, and then threshold it.
But, like a stupid child, i have spent hours in trying to set parameters, use in different combination like: dilate twice, erode three times, blur, and an erode again and so on.
The only result i have obtained is to lost myself: it is not a robust approach at all.
So, for image like this, how can i use those super-popular filters? Or some advice for getting a nice contour of the face and removing the noise?

EDIT:
Those two are two different approach to skin recognition: backprojection and mahalonobis distance, using as training data the output of a simple haar face detection. I'm using both of them because with some illumination setting the first is better, for others the second is better. For my porpouse i can choose the best anytime I'm using my application.

First of all, are you using gray-scale morphology?  If so, then I would threshold the map first, to get a binary image, and then apply binary morphology to fill in gaps and get rid of the noise.

Second, instead of erosion and dilation, try opening and closing, which are higher order morphological operations. Opening is erosion followed by dilation with the same structuring element.  Closing is dilation followed by erosion.

A good thing to try is to do an opening to get rid of the noisy specks, and then do a closing to fill in the gaps.  This will simplify your process to three operations: threshold, open, close.  So you would only need to tune three parameters: the threshold and the two structuring elements.


What are some uses of the overall colour average of an image? As I understand, the overall colour average can be obtained by summing the individual r, g, b values of all pixels within and image and dividing by the number of pixels to produce a new tri-vector.

What are the applications of such techniques?

Well, I can't honestly say that I ever used it, there was a moment when I was considering doing something similar for some quick-and-dirty "shadow" correction.

My friends showed me a flight video, where in the middle of the video there was some dark patches in every frame due to their propeller. To me it looked like a patch in very dark shadow (I never actually got to play with the video so I can't be certain - if the patches actually were the way I suppose or if this would work).

If I didn't want to do any serious correction technique but just spend a bit of time on it for myself, here's what I would do:

calculate the color average of the image $avg_{img}$
find the darker connected areas (the "shadowed" part in each frame) and calculate their average $avg_{dark}$
for all the pixels in the darker area, take their difference from $avg_{dark}$  and add it to $avg_{img}$ and replace their value with this new one
hopefully, result would be a video that's still flawed but somewhat easier to watch than the one with dark propeller-patches flying around in the middle


How to smooth filter response? I have a filter designed in matlab with the function cheby2( N, Rs, Ws, 'stop'). The filter would give nice frequency response with for a given parameter set when the filter order is 2 or 4 (N=4). But if I increase the filter order to say 14 the magnitude plot of filter response is not at all smooth in fact in the stop band it gives very high gain instead of attenuation. Any suggesetion what is going on. What I should look at? Do you need any more information to access the problem. The parameters I am using are:
Rs = 40;
Wc = 0.2;
Wb = 1/128
Ws = Wc +0.5*Wb*[ -1, 1];

You have a very narrow stop band which means that all the poles are crammed in a very small area of the complex plane, close to the unit circle. This can result in severe numerical problems, even for relatively small filter orders, even with floating point arithmetic.

Another important point that you might not realize is that if you design a band pass or a band stop filter using the command cheby2 (and all other similar commands), then N is the order of the prototype low pass filter, and not the order of the resulting band pass or band stop filter. So if you choose N=14 then you've designed a band stop filter of order $28$. Such a filter is very hard, if not impossible to realize because of numerical problems.

Also note that calling cheby2 with [z,p,g] (zeros, poles, and gain) as output argument usually results in higher accuracy than using the polynomial coefficients.

I would actually suggest trying a notch filter instead of such a narrow band band stop filter.


Minimizing sinc spectral lobes for some reasons I'd like to minimize the secondary lobes of the FFT transform of a rectangular signal (shwon in the upper-left side of picture below)

so I've used window function to do that (bottom-left). However, if I want to do it directly on the complex transform of the signal (bottom-right), the values of attenuation of secondary lobe aren't the same as the temporal case.

Having the test code below, I'd like to know what is wrong with it:

pass to the EDIT

     Fs = 1000; L = 2000; t = (0:L-1)*1/Fs;
% the signal
y = zeros(1,L);  y(500:1500) = 1;     % rectangle signal

figure;
subplot(2,2,1);
plot(Fs*t,y); title('Temporal signal'); xlabel('time (milliseconds)')

NFFT = 6*2^nextpow2(L);
Y1 = 2*abs(fft(y,NFFT)); Y1 = Y1/max(Y1);
f = 0:Fs/NFFT:Fs/2;

% the FFT of the signal
subplot(2,2,2);
semilogy(f, Y1(1:NFFT/2+1), '.-'); axis tight; grid on;
title('Signal''s spectrum'); xlabel('Frequency (Hz)'); ylabel('|Y(f)|');

w = window(@blackman, L);
y1 = y.*w';
Y2 = 2*abs(fft(y1,NFFT)); Y2 = Y2/max(Y2);

% the FFT of the windowed signal, in time
subplot(2,2,3);
semilogy(f, Y2(1:NFFT/2+1), '.-'); axis tight; grid on;
title('FFT windowing - temporal product'); xlabel('Frequency (Hz)'); ylabel('|Y(f)|')

% the FFT of the windowed signal, in complex
Y3 = 2*abs(fftshift(conv(Y1, fft(w,NFFT), 'same'))); Y3 = Y3/max(Y3);
subplot(2,2,4);
semilogy(f, Y3(1:NFFT/2+1), '.-');  axis tight; grid on;
title('FFT windowing - spectral convolution'); xlabel('Frequency (Hz)'); ylabel('|Y(f)|')

EDIT:
I've changed the code according to @Matt's remarks, but still there is a problem to the last graph (the noise on the convolutions of the FFTs)

NFFT = 6*2^nextpow2(L);
Y1 = fft(y,NFFT);
YY1 = 2*abs(Y1); YY1 = YY1/max(YY1);
f = 0:Fs/NFFT:Fs/2;

% the FFT of the signal
subplot(2,2,2);
semilogy(f, YY1(1:NFFT/2+1), '.-'); axis tight; grid on;
title('Signal''s spectrum'); xlabel('Frequency (Hz)'); ylabel('|Y(f)|');

w = window(@blackman, L);
y1 = y.*w';
YY2 = 2*abs(fft(y1,NFFT)); YY2 = YY2/max(YY2);

% the FFT of the windowed signal, in time
subplot(2,2,3);
semilogy(f, YY2(1:NFFT/2+1), '.-'); axis tight; grid on;
title('FFT windowing - temporal product'); xlabel('Frequency (Hz)'); ylabel('|Y(f)|')

% the FFT of the windowed signal, in complex
YY3 = 2*abs(cconv(Y1, fft(w,NFFT)', NFFT)); YY3 = YY3/max(YY3);
subplot(2,2,4);
semilogy(f, YY3(1:NFFT/2+1), '.-');  axis tight; grid on;
title('FFT windowing - spectral convolution'); xlabel('Frequency (Hz)'); ylabel('|Y(f)|')'

Multiplication of two time domain sequences corresponds to cyclic (circular) convolution of the DFTs of the two sequences, not to linear convolution as implemented in your code. So if $x[n]$ and $y[n]$ are two length $N$ time domain sequences, and $X[k]$ and $Y[k]$ are their respective DFTs, the following correspondence holds:

$$\text{DFT}_N\left\{x[n]\cdot y[n]\right\}[k]= \frac{1}{N}\sum_{m=0}^{N-1}X[m]Y[(k-m)_N],\quad 0\le k \le N-1$$

where

$$(k-m)_N=k-m+lN$$

with the integer $l$ chosen such that the result is in the range $[0,N-1]$.

Also take a look at this answer to a related question.

EDIT: You also can't use Y1 in the convolution because it is not the DFT but the magnitude of the DFT. You should circularly convolve fft(y) with fft(w).


Unit of Energy Spectral Density The continuous-time Fourier Transform (CTFT) of a signal $x(t)$ (with unit $unit$) is:

$$X(\omega)=\int_{-\infty}^{\infty} x(t)e^{-i\omega t}dt$$

which should be in $unit\cdot sec$ or $\frac{unit}{Hz}$.

The Energy Spectral Density (ESD) of $x(t)$ is defined as:

$$P(\omega)=|X(\omega)|^{2}$$

which should have unit $\frac{unit^2}{Hz^2}$. But that doesn't seem right and certainly doesn't jive with the textbooks.

So what gives?

This question has also been asked here: Power Spectral Density computation and units, but I found the answers unsatisfactory.

If $x(t)$ is an action signal like voltage or current, then the square of $x(t)$ is proportional to instantaneous power and the constant of proportionality depends on what kinda animal $x(t)$ is and what load $x(t)$ is connected to.

So if $\big|x(t)\big|^2$ is proportional to instantaneous power, so also is the integral of $\big|x(t)\big|^2$ over all time proportional to energy (and having the same constant of proportionality).  So you can call this "energy":

$$E_x = \int\limits_{-\infty}^{\infty}\big|x(t)\big|^2 \, dt$$

as long as you remember that it's really proportional to energy and the constant of proportionality depends on what dimension is $x(t)$ and what load that signal is connected to.

For example, if $v(t)$ is a voltage, then $|v(t)|^2$ is proportional to instantaneous power and if $v(t)$ is connected to a resistor having resistance $R$, then instantaneous power is

$$ p(t) = \frac{1}{R}\big|v(t)\big|^2 $$

and the total energy is

$$E_v=\int\limits_{-\infty}^{\infty}p(t) \, dt$$

or

$$E_v=\frac{1}{R}\int\limits_{-\infty}^{\infty}\big|v(t)\big|^2 \, dt$$

Clearly the constant of proportionality is $\frac{1}{R}$.

If $x(t)$ is a dimensionless value in a DSP or computer, then you have to scale it with the reference voltage $V_\mathrm{ref}$ of the D/A converter to make this number a voltage that delivers power to a load.  Then the constant of proportionality becomes $\frac{V^2_\mathrm{ref}}{R}$ and the energy really is (with correct scaling):

$$E_x = \frac{V^2_\mathrm{ref}}{R} \int\limits_{-\infty}^{\infty}\big|x(t)\big|^2 \, dt$$


Autocorrelation function of a discrete signal For a certain experiment of which I have data, I want to calculate the autocorrelation of the measurements. For this I used the expression from Wikipedia, which reads

Here, $n$ is the dataset's size, $0 < k < n$, $\sigma$ is the standard deviation and $\mu$ is the mean.

Now, I've read the descriptions on the page, but I don't fully understand all of the terms in there. I understand why we divide by $(n-k$), because the number of datapoints used constantly changes (you have a finite signal, after all). However, I don't see why we divide by the variance. Moreover, I also don't know why we subtract the mean. Doesn't this lead to problems, if for example for n = 50 and k = 49, we want to multiply $(x_1-\mu)(x_{50}-\mu)$, and $x_1 > \mu$ but $x_{50} < \mu$? I don't really see what a negative correlation would entail.

For my data, which is only 0's and 1's, this is even more problematic. Using the above formula the correlations quickly decay to zero, but at the end they violently start to spike, because the contribution from the mean becomes so big, relatively. In my data the mean is 0.5 and the variance is 0.25, so for example, my final point $R(49)$ = $\frac{(0-0.5)(1-0.5)}{0.25}$ = -1

But -1 is huge, especially when 'previously' the correlations had died out to 0. How does one deal with this?

An illustration of this is

First: you should really read some basic theory of autocorrelation. For example.

  I also don't know why we subtract the mean

First, to subtract the mean is the usual and right thing to do – it's so standard that often it is straightly assumed that the signal has zero mean. Recall that the covariance of two variables is $E[ (x-\mu_x)(y-\mu_y)]$ – of which the variance is a special case. Here we are assuming that the process is stationary, so the mean is constant.

  I don't really see what a negative correlation would entail.

The correlation (as the covariance) can be positive, zero or negative. A negative correlation says that the values tend to have opposite values, that is: when X has a "big value" (relative to the mean), then it's more probable that Y has a "small value" (relative to the mean). If the correlation is positive, then's the reverse.

  However, I don't see why we divide by the variance.

You don't have to. You divide by it if you want a normalized correlation (rather: a correlation coefficient), that will be in the $[-1, 1]$ range.

  I understand why we divide by $(n−k)$...

Actually you have two alternatives here: if you divide by $(n-k)$ you get an unbiased estimator (good), but, you get large variance (the "noise" you see in the graph) for large values of $k$ (bad; however, one is often interested only in the correlation for small values of $k$). Alternatively, one can divide by $n$; this produce a biased estimator but better behaved (lower variance) for large $k$. This is also explained here or in any textbook about correlation estimators.


Frequency response of system with feedback and feedforward co-efficients? For seismological application, I am playing with an elliptical filter, and once I have both sets of feedback and feedforward co-efficients, how do I make it so that I can plot the frequency response of the system? Can I also use the same technique to make the frequency response for a bessel?

I know that I can take the FFT for an FIR filter and I can get its transfer function, but for a non-FIR, is there a way also using the FFT?

Yes, you can always use the FFT to compute the frequency response of a discrete-time system at discrete frequencies. For an IIR filter you get the following transfer function:

$$H(z)=\frac{\sum_{m=0}^{M}b_mz^{-m}}{\sum_{n=0}^{N}a_nz^{-n}}\tag{1}$$

Usually, $a_0$ is normalized to 1. In order to evaluate $H(z)$ on the unit circle $z=e^{j\theta}$ you can compute the FFTs of the numerator and of the denominator of (1) separately and then divide to get the frequency response of the filter.

If you use Matlab or Octave, you can simply use the command $\tt{freqz}$.


estimation of ARMA(1,1) with MA parameter greater than 1 I am working with the following simple ARMA(1,1) model:
$$
z_{t+1} = \phi z_{t} + \theta\varepsilon_{t} + \varepsilon_{t+1}
$$
In my case $\theta$ depends on some other parameters, and, therefore, I know that $\theta > 1$.

Let us generate 1000 ARMA(1,1) samples with $\phi = 0.95$, $\theta = 0.5$,
$\sigma = 0.08$  and estimate the parameters.

import numpy as np
import statsmodels.api as sm
from statsmodels.tsa.arima_process import ArmaProcess

process = ArmaProcess(np.r_[1, -0.95], np.r_[1, 0.5])
y = process.generate_sample(1000, scale=0.08)
model = sm.tsa.ARMA(y, (1, 1)).fit(trend='nc', disp=0)

The result is [0.96069232 0.51912881 0.080], which is quite close to the true values.

Next, lets us try parameters $\phi = 0.95$, $\theta = 2$,
$\sigma = 0.08$, i.e. change the following row in the code above:

process = ArmaProcess(np.r_[1, -0.95], np.r_[1, 2])

The result is [0.95668055 0.49700497 0.156] and there is no errors and warnings!

By default, in order to have distinguishable auto-covariance function we assume $|\theta| < 1$. Therefore, in the result the estimator is close to $\frac{1}{\theta}$. Though, the standard deviation of the variance is far from the true value.

Is there a way to estimate the case when $\theta > 1$? From what I can see, all packages use the fact that $|\theta| < 1$.

I suspect that the reason is that many ARMA fitting algorithms require that the ARMA system be stably invertible.  In your example, the MA process is not minimum phase as the zero will be outside the unit circle. That makes the inverted system unstable (the pole will be outside the unit circle).

I don't know the details of the sm.tsa.ARMA(y, (1, 1)).fit algorithm, but I suspect that it will also be so restricted.

See, for example, this lecture by Florian Pelgrin, which contains the slide below.

In particular, note the second condition.

This means that I suspect it will be harder to find an algorithm that does the estimation for the case you're  interested in.


How do I find the power of a specific frequency at a specific sample of a discrete time varying signal? EDIT:

First off let me explain that I understand that frequency makes no sense for a single sample. What I am actually talking about is the power spectrum within a short time window around a specific sample.

I have a time varying discrete signal: $\ s(t_i), i=0...199 $ with sampling dt = 10 ms.

I want to find the power of a specific frequency: 15 Hz at time sample j = 100.

This is how I would solve this:

I understand that I must extract a subset of the signal in a window around j = 100 and multiply it with a window function to avoid "ringing".

I therefore extract the subset j=93..107 and multiply it with a -/+ 2 standard deviations gaussian kernel.

Any advice on the size of the window?

Next I zero-pad the subset to 512 samples in order to achieve some "spectral interpolation"
and take the FFT.

Finally I locate the complex number corresponding to f=15 Hz and take the absolute value.

Actually I would repeat this procedure for j = 7...192 to find the power of f=15 Hz at each sample.

However reading about the Short Time Fourier Transform got me confused:

they seem to break the data to be transformed into overlapping chunks,
window the data within each chunk, take the FFT and then somehow adding the results together!?

Is this just a matter of performance?
As a first solution I would prefer a simple but slow approach (such as mine).

I realize this doesn't answer the detail of what you're asking, but I believe it answers the question in the title:

See Rick Lyons's writeup here for using the Goertzel algorithm to get the signal level at a specific frequency.

You can then follow up the equations:

with:

$$
\hat{p}(n) = | y(n) |^2
$$

and use $\hat{p}(n)$ as your "instantaneous power" at frequency $\omega =2\pi m/N$ --- which is sort of like  Lyons's eq. 13.83:

(though he's evaluating this one at the "end" of his DFT length).


why this frequency response is added to the power of e? From where did this ${\angle}H{(\Omega)}$ came into power of $e$ ?

$H(\Omega)$ is a complex number. $Ae^{j(\Omega t + \phi)}$ is another complex number, expressed in polar format by its magnitude $A$ and phase $\Omega t + \phi$. When you multiply complex numbers in this format, you simply multiply their magnitudes and sum their phases. Thus:

$$
H(\Omega)Ae^{j(\Omega t + \phi)} = |H(\Omega)| Ae^{j(\Omega t + \phi \angle H(\Omega))}
$$


How to custom optimize cuFFT for a mini batch of multi-channel images? I am reading a paper, which has the following paragraph.

  Current GPU implementations of the FFT such as cuFFT are designed to
  parallelize over individual transforms. This can be useful for
  computing a limited number of transforms on large inputs, but is not
  suitable for our task since we are performing many FFTs over
  relatively small inputs. Therefore, we developed a custom CUDA
  implementation of the Cooley-Tukey FFT algorithm which enabled us
  to parallelize over feature maps, minibatches and within each 2-D
  transform. Note that 2-D FFTs lend themselves naturally to
  parallelization since they can be decomposed into two sets of 1-D FFTs
  (one over rows and the other over columns), and each set can be done
  in parallel.

Unfortunately the paper doesn't bring further detail on this. As I could not get a feedback from the authors, I am asking the community,

How this implementation differs from the standard cuFFT implementation of the FFT?

The only thing that comes to my mind is that the standard cuFFT for images, when applied to a mini-batch, performs individual 2D FFT. Then, the only re-computation that could possible be avoided when considering a mini-batch of multi-channel images would be the re-calculation of the Twiddle factor $W_{N}^{k}$. Would there be another way to speed up a 2D FFT for this scenario?

I don't know if the paper is right, but what it claims is pretty clearly stated there:

cuFFT uses many GPU processors in parallel to calculate one large FFT fast. It splits a complex problem in many smaller problems and solves them in parallel, and then combines the results to the solution to the complex problem.
Their implementation goes and calculates many small transforms in parallel. It takes many small problems and solves them in parallel.

To be honest, that in isolation is not that cool, considering all they had to do is take a single-threaded FFT and tell CUDA to execute it on different data in parallel.

That trivial data parallelism is exactly what CUDA does for you. There were implementations of that even before Nvidia had a stable cuFFT out there. In fact, multiple parallel FFTs done on a GPU have been around longer than CUDA – it's the obvious algorithm to implement on a manycore thing like a GPU.

I doubt the authors are fully right in their claim that cuFFT can't calculate FFTs in parallel; cuFFT especially has a function cufftPlanMany which is used to calculate many FFTs at once.

To cite the cuFFT documentation:

  where batch denotes the number of transforms that will be executed in parallel,

so, either that documentation is wrong, or the authors simply implemented something that was already built-in to cuFFT but which they ignored (maybe for good reason). Maybe that functionality wasn't there in 2013 (but I'm not willing to bet on that).


Initial rest condition for the linear constant-coefficient differential equations Suppose that system has the input/output relation as follows

$$\sum_{k=0}^{N}a_k \frac{d^ky(t)}{dt^k} = \sum_{k=0}^{M}b_k \frac{d^kx(t)}{dt^k}$$
Where $a_k, b_k \in \mathbb{R}$.

Obviously we need additional conditions to determine completely the input-output relationship for the system and those auxiliary conditions specify whether the system is LTI and causal or not.
The definition of initial rest condition for a system in general is that if $x(t) = 0$ for $t\lt t_0$ then $y(t) = 0$ for $t\lt t_0$.
First of all is it true that the initial rest condition for the all systems implies $y(t_0) = 0$? Also I have seen expressions like $y({t_0}^+)$ and $y({t_0}^-)$ but I don't know how they relate to $y(t_0)$. Another question is that whether the initial rest condition is necessary and sufficient for the mentioned system to be causal and LTI? What's the rigorous proof, then? I searched on the Internet but didn't found a proof or a clear explanation.

The notion of initial rest and causality follows only for linear systems (time invariant not required).

For proof of Initial rest being sufficient condition for an LTI system to be causal. Please refer to Problem 1.44 of Allen Oppenheim's book on Signal and Systems. Try and solve it yourself (mostly a homework problem, there is also a solution manual just to let you know in case you are stuck), it doesn't require anything extra other than convolution sum.

For $y(t_o^+)$ and $y(t_o^-)$ relating to $y(t_o)$ post the quote you are looking for clarifications on, it depends on multiple context, generally that would be specifying an initial condition of  a one sided Fourier transform or one sided sequence.


Mathematical Proof of Complementary Filter？ Having found some unofficial sources on Complementary Filter (Thousand Thoughts Sensor Fusion and The Balance Filter by Shane Colton), I wish to work out its rigorous mathematical proof.

The Complementary Filter,
$$y=\alpha \times y+(1-\alpha) \times x$$
where $\alpha$ is the filter parameter, usually chosen to be ~0.98, is named as such, because effectively the filter highpasses $y$ and lowpasses $x$.

With this setup ($\alpha=0.98$), man claims that the filter can filter out the low-frequency part of $y$ and meanwhile the high-frequency part of $x$.

What is the mathematical proof here?

There is no rigorous proof so far available. All the sources simply take it as granted and validates it with experiment data, which is unacceptable in scientific validation.

This is a first order, low pass infinite impulse response filter (IIR).  It is some times called the exponentially weighted moving average.  You can find the derivation on Wikipedia. Orfanidis' Introduction to Signal Processing gives a nice description as well.

Don't think of it as high passing y and low passing x, this will lead to confusion.  It is weighting y more heavily than x, the consequence is short time variations of x have a small effect on y, therefore removing higher frequencies.

Also looking at your other questions, it appears that you think that this is a high pass filter as well.  This is not the case.  α is valid from 0 to 1, so you can see that you will get a lot of low pass filtering (High time constant, limited by your digital precision) or no filtering where y = x ( In the case of x = 0).  Read up on high pass filters if you need one, wikipedia gives the difference equation (you have provided the low pass one) for a single order high pass filter.


Calculating the first derivative of an image using DFT I need to calculate the first derivative of a greyscale image (a 2D array) using a DFT function I built (which works). Unfortunately, the results don't seem to be correct.

In the fourier domain, the derivative d/dx is given as F(u,v)* 2*pi*i/N * u, where u is the x-axis transformed, N is the size of one of the matrix's dimensions, v being the other one.

Attached is the code. What bothers me is that I'm not getting the same results as I would differentiating by convolution with (1,-1) or (1,-1) as a column vector.

def derivative(fourier_signal):
"""
Derivative in fourier domain is multiplying by u or v, and 2pi*i/N
:param fourier_signal:
:return:
"""
N = np.shape(fourier_signal)[ZERO]
M = np.shape(fourier_signal)[ONE]
u = np.arange(N)
v = np.arange(M)
du = fourier_signal * (u*TWO_PI*1j)/N
dv = fourier_signal * (v*TWO_PI*1j)/M
return du, dv

def fourier_der(im):
    # Calculate DFT2
    dft_image = DFT2(im)
    # Function that Multiply by rows by u, columns by y
    du, dv = derivative(dft_image)
    shifted_du, shifted_dv = np.fft.fftshift(du), np.fft.fftshift(dv)
    dx, dy = IDFT2(shifted_du), IDFT2(shifted_dv)

I'm not looking for easy answers on how to do it, but rather a direction to as to why my output is incorrect.

They are not the same. Using 1D notation,

the discrete-time (backward) first difference is $x[n] - x[n-1]$ whose frequency domain DTFT equivalent is

$$ x[n]-x[n-1] \leftrightarrow X(e^{j\omega}) - e^{-j \omega} X(e^{j\omega}) =X(e^{j\omega})(1- e^{-j \omega})  $$

which becomes
$$ x[n]-x[n-1] \longleftrightarrow X[k](1 - e^{-j \frac{2 \pi}{N} k})$$ using the DFT to implement it.

The FIR impulse response of the discrete-time system that implements the first difference is therefore,

$$ h[n] = \delta[n] - \delta[n-1]$$

The first derivative of a continuous-variable function $x(t)$ is $x'(
t)$ and in CTFT domain it becomes :
$$ x'(t) \longleftrightarrow j\Omega X(\Omega) $$

where the analog system frequency response is

$$H_d(\Omega) = j \Omega $$

which is not implementable in digital form, but a bandlimited approximation to it is attained under a sampling period of $T$ that yields an equivalent discrete-time frequency response of a (bandlimited) differentiator as

$$ H_d(e^{j \omega}) = j \frac{\omega}{T} ~~~, ~~~~ |\omega| < \pi$$

The associated discrete-time (IIR) impulse response is
$$ h_d[n] = \text{IDTFT} \{ H_d(e^{j \omega}) \} = \text{IDTFT} \{ j \frac{\omega}{T}  \} $$

Practically you would truncate and window $h_d[n]$ before using.

So they are not the same.


Differentials - Differences: Non causality in the system I'm still learning DSP and referring to Oppenheim video lectures.

In that lectures, differential difference equation is obtained for IIR filter design, in Lecture 14.

$$\mathcal{L}[\frac{\mathrm d}{\mathrm dx}y_a(t)] = sY_a(S)$$
$$\mathcal{Z}[\frac{y[n+1] - y[n]}{T}] = \frac{z-1}{T}Y(z)$$

And using these two equations, difference equation $z-1=sT$ is obtained.

To calculate the z-transform of the differential, $y(n+1)$ is used. By using that term, are we making our system non-causal..? (As our equation depends on future values)

Thank you for your time..

Replacing the derivative by forward differences is just a way to define a transformation from the $s$-plane to the $z$-plane, and it has nothing to do with the causality of the resulting discrete-time filter. As an exercise you could try to do the same with backward differences.

Note that this method for transforming a continuous-time filter to a discrete-time filter is useless in practice because stability is generally not preserved, and the transformed frequency response may be drastically different from the original response.


Topological skeleton matching techniques What are the most common and best topological skeleton matching techniques?

Suppose I have skeletonized object and corresponding 2D model of this object. What are the next steps to match pattern? What are skeleton descriptors?

If someone could give code examples, reference to some open source projects using skeletonization and matching I would be grateful.

I want a method that will work also when skeleton is lack of some part of it due to bad object detection, but any suggestions are appreciated.

I think you should look/google in the direction of graph matching.

This is because you can represent your skeleton as a graph by creating nodes at endpoints and crossings, and creating edges for the connected endpoints and crossings. There is much more research done in the direction of graph matching than just topological skeleton matching. However for the latter I found this paper, which you might find interesting.

You can also look through the references of the Wikipedia page of Topological skeletonization, there are also some references about matching skeletons.


Detecting a frequency pattern in the frequency domain I am looking to detect a combination of frequencies in the frequency domain. For example, here is my spectrum:

And I'm trying to detect the pilot tone. I have the FFT of the pilot tone:

which zoomed in looks like this:

I tried doing a cross-correlation in matlab:

trainingfft = fft(trainingfreq);
datafft = fft(idat);
Hxcorr = dsp.Crosscorrelator;
CorrResult = step(Hxcorr,datafft,trainingfft);

But ended up with this:

I really am looking to detect the presence of the three peaks anywhere in the spectrum. I think the cross correlation isn't working because of the presence of the noise floor, so I was thinking about using some kind of threshold checker to clean up the noise but really am looking for an elegant solution.

Is there an elegant way of doing this?

Try

xcorr(trainingfft,datafft,'option')

to specified a normalization for the cross-correlation. See the file page for more details of the option.


Deconvolution of Images - How To? I'd like to do a deconvolution of image. For example for convolution I'm using a $3\times 3$ mask with all elements $= 1$:

$$\begin{bmatrix}1 & 1 & 1 \\ 1 & 1 & 1 \\ 1 & 1 & 1\end{bmatrix}$$

For a pixel $X$ as follows:

$$\begin{bmatrix}X_1 & X_2 & X_3 \\ X_4 & \ \mathbf X & X_5 \\ X_6 & X_7 & X_8\end{bmatrix}$$

which is convolved by a kernel given by:

$$\begin{bmatrix}a & b & c \\ d & e & f \\ g & h & i\end{bmatrix}$$

I believe $X$ should be determined as follows

$$X=\frac{aX_1+bX_2+cX_3+dX_4+eX+fX_5+gX_6+hX_7+iX_8}{a+b+c+d+e+f+g+h+i}$$

This is a basic convolution. Now I'd like to do a deconvolution of an image which was convolved with above mask. Does anybody know how to count that?

I will try to address @Hilmar answer and answer the general question.

Hilmar Answer

  Deconvolution is in the general case not possible, so it needs to be
  approximated with application specific constraints and requirements.

  Let's look at a simple 1-dimensional example that illustrates the
  problem. Assume you have an impulse response like your kernel, i.e. h
  = [1 1 1]. Then let's look at the output of a signal (which you can think of a line of pixels), x = [2 -1 -1 2 -1 -1 2 -1 -1 2 -1 -1 2 -1
  -1]. The result of the convolution is

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

The case above is simply wrong. Specifically if one looks on my comment:

  The result being zero doesn't mean data is lost. It has the same information as any other number.

Let's do the following example:

% Input Signal
vX = [2; -1; -1; 2; -1; -1; 2; -1; -1; 2; -1; -1; 2; -1; -1];
% Convolution Kernel
vH = [1; 1; 1];

% Convolution Result
vC = conv(vX, vH, 'full');

% Convolution Operator - Matrix Form
mH = CreateConvMtx1D(vH, length(vX), 1);

% Deconvolution
vXR = mH \ vC;

Where CreateConvMtx1D() is given in my GitHub Code.

The result is given by:

vXR.'

ans =

    2.0000   -1.0000   -1.0000    2.0000   -1.0000   -1.0000    2.0000   -1.0000   -1.0000    2.0000   -1.0000   -1.0000    2.0000   -1.0000   -1.0000

Which is a perfect reconstruction of the input signal.

General Case

Indeed the properties of Deconvolution in 1D and 2D are similar.
The ability to reconstruct the signal are basically a function of 2 parameters:

The SNR of the data (Added noise).
The Condition Number of the Convolution Operator of the problem.

Let us assume there is no added noise in the problem.
Then we're only limited by the Condition Number.
In theory, for any bounded Condition Number one would be able to perfectly reconstruct the signal. Yet since we use finite and quantized representation of numbers (Let's say Float64) we are limited with the condition number we're able to tackle.

For more details see - 1D Deconvolution with Gaussian Kernel (MATLAB).


How do I estimate the derivative of the current position, when I have only values from past to present? If I have a discrete real-time signal $x[n]$, with its latest value $x[i]$ and all its past values $x[i-t]$, how can I estimate the derivative at $x[i]$?

Is the derivative between x[n] and x[n-1] useful enough, hence x'[n] = x[n] - x[n-1] ? Or do you need a better estimate?


How to plot ROC curve for weka multilayerperceptron? I am working on an image classification project I use multilayerperceptron,how do I plot ROC curve.I use weka's perceptron.

From University of Waikato WEKA wiki:

  The class used for calculating the ROC and also the AUC (= area under the curve) is weka.classifiers.evaluation.ThresholdCurve.


Why can we have non-integer frequency bins in FFT? I am studying DFT/FFT and I'm very confused about one thing. I read online that the frequencies we can sample with DFT must be integer (Why does the frequency in the DFT have to be an integer?). Later I found out that the frequency bins in frequency domain in FFT can be non-integer numbers, since they are calculated like this:

delta f = sampling frequency / number of samples

Why can the frequency bins be the non-integer number here? Can we also sample non-integer frequencies using DFT/FFT?

Thank you for your help.

Frequencies in DFT, $\omega = \frac{2\pi}{N}k, \ k = 0,1,2,...,N-1 $, only depends on the Length $N$ of DFT, and nothing else.

$$X[k] = \sum^{N-1}_{n=0}x[n]e^{-j\frac{2\pi}{N} nk}, \qquad k = 0,1,2,...,N-1$$

It is very important to understand this expression as the projection of time-domain finite length sequence $x[n], n=0,1,2,3,...,N-1$, of length $N$, onto orthogonal Discrete Fourier Basis Vectors $W_k^{N} = e^{j\frac{2\pi}{N}nk}, k=0,1,2,...,N-1$. And, the DFT Coefficients $X[k]$ are the projection coefficients actually on the $k^{th}$  basis vector. So, with this idea in mind, we can say $x[n]$ is actually a vector $\vec{x}$ in an $N$-Dimensional vector space. And, this $N$-D vector space has a set on $N$ Fourier Basis which reveals the frequency information of that $\vec{x}$. We know that projection of a vector onto another depends on dot-product (inner-product), $\langle\vec{x},\vec{y} \rangle$. The projection coefficient of a vector $\vec{x}$ onto a vector $\vec{y}$ :

$$P_{\vec{x}} = \frac{\langle\vec{x}, \vec{y}\rangle}{\langle\vec{y}, \vec{y}\rangle} = \frac{\sum^{N-1}_{n=0}x[n]y^{*}[n]}{\sum^{N-1}_{n=0}y[n]y^{*}[n]}$$
And, if we are projecting a vector $\vec{x}$ onto a unit length vector $\hat{y}$, then the Projection coefficient is exactly the dot-product, because, the denominator of the above expression will become 1.

$$P_{\vec{x}} = \langle\vec{x}, \hat{y}\rangle$$
That is exactly what the DFT Coefficient is. Projection Coefficient of time-domain sequence $\vec{x}$ onto $\vec{W_k^N}$
$$X[k] = \langle\vec{x}, \vec{W_{k}^{N}}\rangle = \sum^{N-1}_{n=0}x[n]\,W_k^{N}[n] = \sum^{N-1}_{n=0}x[n]e^{-j\frac{2\pi}{N} nk}$$

Now, notice that the $N$ Fourier Basis vectors are each at a digital frequency $\omega = \frac{2\pi}{N}k, k=0,1,2,...,N-1$. So, the DFT/FFT frequency bins just depend on the DFT Length and nothing else. DFT by itself is a perfectly fine operation on an $N$-dimensional vector. It is a linear transformation, just a change of basis from canonical basis set to Fourier Basis set, so as to reveal frequency information in that signal $\vec{x}$.

It reveals frequency information because projection(inner-product) of one vector onto another actually is a measure of similarity between the two vectors. So, each of the $X[k]$ is telling us how much similarity is between the time-domain sequence $\vec{x}$ and $k^{th}$ Fourier Basis, and, $k^{th}$ Fourier Basis is a sequence of digital frequency $\omega = \frac{2\pi}{N}k$.

As I said, DFT by itself is perfectly fine operation on any vector of finite length. You do not have to link these digital frequencies with sampling frequency. But if you want to, then, it is enough to know that when we sample a signal at sampling frequency $f_s$, then the maximum frequency that can be represented in the sampled signal is $f_s/2$. And the maximum digital frequency is $\pi$. Hence, $f_s/2$ maps to $\pi$. And therefore the digital frequency of $k^{th}$ Fourier Basis, $\omega = \frac{2\pi}{N}k$ maps to $f = \frac{f_s}{N}k$.


compression ratio of pixels block I am new to image compression topic and really appreciate it if someone could help me understand how compression rati works in this case. Assuming I have JPEG lossless coder that encodes the following image:

2 3 4 3
1 A B C

I have used predictive coding and Huffman to obtain the binary representations (codes) for pixels A, B and C:

Which are:

A=0 -> 0
B=3 -> 1101
C=0 -> 0

How can I calculate the compression ratio achieved for coding of the pixels A, B and C (jointly) assuming that the original image has pixel values ranging from 1 to 8?

Compression ratio is $$ \frac{N_{uncomp}^b}{N_{comp}^b}$$
where ${N_{comp}^b}$ is the total number bits requried to represent those $3$ pixels $A$,$B$, and $C$ while ${N_{uncomp}^b}$ is the total number of bits required when they are not compressed.

Based on your claims, assuming that, those three pixels had $8$ levels with fixed length coding (FLC) of $3$ bits per pixel then; ${N_{uncomp}^b} = 3 \times 3 = 9$ bits in total. And again based on your variable length code (VLC) of the lossless coding scheme which uses $1$ bit for $A$ and $C$ and $4$ bits for $B$ you have ${N_{comp}^b} = 1 + 4 + 1 = 6$ bits in the compressed case, hence the ratio is $$ \frac{N_{uncomp}^b}{N_{comp}^b} = \frac{9}{6} = \frac{3}{2} $$

which means a $1.5 \times$ reduction in the number of bits required to represent the same source compared to its uncompressed form.

Note that you can also use the reciprocal of that ratio to indicate a similar metric, which would be $2/3$ in this case indicating the size that compressed form will have compared to the original.