Keeping track of units for spectral energy of discrete signals The energy of a continuous signal $x(t)$ recorded between $t=0$ and $t=T$ (namely $x(t \notin [0, T]) = 0$) is defined as

$$E = \int_0^T |x(t)|^2\, dt.$$

In most signal processing texts I met, the formula for the energy of the same signal sampled at $N$ regular steps separated by $\Delta t$ is given as

$$E_N = \sum_{n=0}^{N-1} |x_n|^2,$$

which implies an implicit normalisation of the signal.
As I would like to keep track of the physical units, I am inclined to use

$$E_N = \sum_{n=0}^{N-1} |x_n|^2 \Delta t.$$

Now, Parseval's identity gives, in the continuous case,

$$E = \int_0^T |x(t)|^2\, dt = \int_{-\infty}^\infty |\hat{x}(f)|^2\, df,$$

where $\hat{x}$ is the Fourier transform of $x$.
This leads to define the /energy density/ as

$$e(f) = |\hat{x}(t)|^2.$$

I would translate this in the discrete case as

$$E_N = \sum_{n=0}^{N-1} |x_n|^2 \Delta t = \Delta t^2 \sum_{k=0}^{N-1} a_k^2 \Delta f,$$

where $\Delta f = 1/T$, and the $(a_k)$ are the Fourier amplitudes and such that

$$\forall \; k \in \{0, N-1\}, \, \sum_{n=0}^{N-1} x_n e^{-2i\pi k \frac{n t}{T}}
                                      = a_k e^{i \phi_k},$$

Finally, I would define the discrete energy density as

$$e_N(f) = \Delta t^2 |a_k|^2 \text{ for } k \Delta f \leq f < (k+1) \Delta f.$$

Does this sounds sane to practitioners of the field?

First it should be noted that the nomenclature using power and energy as properties of signals is somewhat misleading. For physical energy to come out of a time integral over a squared signal, the signal would need to be of dimension $\sqrt{[\mathrm{energy}]/[\mathrm{time}]}$, which is not the case for practically all signals that you encounter.

In fact, if you have signals that fully describe the state of a physical system over time, and the energy of the system can there be expressed using just that signal, the result is usually very different and most importantly, it is a property of each time instant and not the whole signal. If you have a signal that describes a non-temporal aspect of a physical signal, like the momentary state of a string, the energy of the system will in general require the time derivative of that state.

Therefore, the terms energy and power are a bit of stretch and should be used carefully. I personally prefer speaking of a squared magnitude and a squared magnitude density.

With this preamble, the use of physical units like time and frequency in the definition of power and energy of a signal become somewhat relative. The choice of scaling and therefore of units is mostly arbitrary and depends on what you want to express and which transformations you want to keep your measures invariant.

For example, if you have a choice of sampling your physical system at different rates and you want the energy measure to be invariant under such rate changes, then you need to first understand how the amplitudes behave under a rate change. The normal behaviour of an amplitude is to be be unaffected by a change of sampling rate. Some acquisition methods may however rely on integration to record the samples and your sample values may therefore scale with the sampling rate. So your energy sum may contain different powers of the time step in order to keep the energy invariant under sampling rate changes.

Keeping track of all this for different kinds of signal behaviours quickly become tricky. That is why I made a habit of defining a discrete signal not just by giving a sequence of numbers, but also by providing a homomorphism to an infinite dimensional function space that allows me to work out all the meaningful measures and properties in a well defined framework with properly defined integral transforms and continuous variables.

So, a discrete signal is a tuple $\{ s\in \mathbb{S},\, D: \mathbb{C}\to\mathbb{S},\, \bar{D}:\mathbb{S}\to\mathbb{C} \}$ where $\mathbb{S}$ is the space the discrete signal lives in, $\mathbb{C}$ the space of the continuous time signals, and $s$ a representation of the discrete time signal in $\mathbb{S}$. Note that $D\circ\bar{D}=\mathrm{id}_\mathbb{S}$ but not necessarily $\bar{D}\circ D=\mathrm{id}_\mathbb{C}$.

You give these two maps $D$ and $\bar{D}$ all the invariance properties you would like to see. Afterwards you can pull back all definitions from $\mathbb{C}$ to $\mathbb{S}$ by applying the two domain changing maps.

This method of working with discrete signals will give you a well defined way of talking about properties define for continuous signals and allow you to get your scaling correct under sample rate changes or changes of other parameters relevant for your discrete signal representation.

I will not try to replicate your derivation using this approach, because it is not clear to me which invariance and signal properties you want to preserve and what exactly your goal is. Like I said, it depends on the context what is right. I still hope this answer helps your understanding and answers the slightly more fundamental question of how to deal with units, dimensions and scaling of discrete signals compared to continuous signals.


Can a pole at the origin be considered a dominant pole? I am studying control syetems, and I have studied that the dominant poles of a system are the ones closer to the imaginary axis.

Now, consider a system with a pole at the origin. Can this pole at the origin be considered a dominant pole?

You could get away with defining it as such (or extending the concept to the pole with the most positive real part).

But typically the concept of a dominant pole is applied to systems that are truly stable, and thus actually settle out.  A system with singular poles on the imaginary axis (and no poles in the right-half plane) is considered "metastable", because with the right input it has a bounded response, but with the wrong input its response is unbounded.


Fourier Transform Form: two sin components & a phase shift & a magnitude for only one term This is an example from my text book of a continuous signal:
$$x_{in}(t)=\sin \left( 2\pi \cdot 1000 \cdot t\right) + 0.5\sin \left( 2\pi \cdot 2000 \cdot t + \dfrac{3\pi}{4} \right) $$
So to perform a fourier transform on this signal, how to do that, isn't it a bit funny, since it has two sine components. Shouldn't complex numbers have a sine term and a cosine term? And it's got a scalar term applied to only one component, don't those usually apply across both terms of a complex component? And it's phase shifted, what to do about that?

Fourier Transform is a linear one, so you can make use of superposition principle:

$$ \mathscr{F} [\alpha x(t) + \beta y(t)] = \alpha \mathscr{F}[x(t)] + \beta \mathscr{F}[y(t)] $$

So for the first component $$x(t) = \sin \left( 2\pi \cdot 1000 \cdot t\right)$$

by definition:

$$\mathscr{F}\left[\sin(2\pi f_0 t + \phi) \right] = \dfrac{i}{2} \left[ e^{-i \phi}\delta(f+f_0) - e^{i \phi}\delta(f-f_0) \right] $$

you get:

$$ \mathscr{F}[x(t)]=\dfrac{i}{2} \left[ \delta(f+1000) - \delta(f-1000) \right]  $$

Second component is a sinusoid with shifted phase, so the complex exponent represents that:

$$y(t) = \dfrac{1}{2} \sin \left( 2\pi \cdot 2000 \cdot t + \dfrac{3\pi}{4} \right)$$

has following Fourier Transform:

$$\mathscr{F}[y(t)] = \dfrac{1}{2}\dfrac{i}{2} \left[ e^{\dfrac{-3\pi i}{4}}\delta(f+2000) - e^{\dfrac{3\pi i}{4}}\delta(f-2000) \right] $$

By summing both results you get the Fourier Transform of your signal.


spectral subtraction for stationary noise my recording was suffering from stationary background noise. I used spectral subtraction to remove it. the noise is reduced in loudness but not eliminated completely.
before:
http://vocaroo.com/i/s18FtoUwSah4
after:
http://vocaroo.com/i/s12U4feLs6iu

is spectral subtraction working?What can I do to further improve my results.

In almost all practical applications of spectral subtraction it's unrealistic to expect that the noise can be completely eliminated. Even if this were possible, the price to be paid would be severe distortion of the desired signal.

What you can try is tweak the parameters. The most obvious parameter to tune is the (over)-subtraction factor, which determines how much of the estimated noise is subtracted from the noisy spectrum. It is not uncommon to have over-subtraction factors of around $2.0$ or $2.5$. Of course you have to make sure that the resulting estimated speech (magnitude or power) spectrum remains non-negative.


Detecting the original bandwidth of upsampled recordings I have a set of recordings sampled at 16khz, but some of them are upsampled from 11khz and some from 8khz. Quite some of the recordings are noisy. What is the best way to classify the recordings according to the original bandwidth?

I know some people are using GMM classifiers for detection. A paper on the subject is DETECTING BANDLIMITED AUDIO IN BROADCAST TELEVISION SHOWS

I've come to the following algorithm: calculate energy for 40 frequency bins for the full range from 0 to 8khz for 10ms frames, then normalize by maximum energy, then calculate maximum of normalized energy over sequence of frames. Bandwidth is detected by comparing maximum normalized energy in high-order bins compared to maximum normalized energy in middle bins.

Any better ideas?

A recording originally at 8kHz and digitally upsampled to 16kHz will have almost no energy in the 4-8kHz range (whatever is here is due to imperfections in the filters used for the upsampling process). I would just use a 4kHz and 5.5kHz high pass; and use a threshold on the signal energy at the output of these filters.

... Unless your recordings are affected by noise, distortion etc. after they have been upsampled to 16kHz. Is that the case? If so, please update your question and tell the whole story about your signals!


Remove spurious peaks of Hilbert transform I want to eliminate spurious peaks of Hilbert transform for finding glottal closure in the linear prediction (LP) residual. I have the following 4 steps:

Down-sample.
Hilbert Transform.
Identify Peaks in Hilbert Transform.
Consider this hypothesis that time gap between two successive glottal closure instance is not likely to vary in the adjacent pitch period.

Case 3 would help to find peaks such as candidate and spurious peaks. However case 4 improve our results for best peaks candidate for glottal closure.

I have written the following MATLAB code:

e_downsample = downsample(e,2);
DFT=fft(e);
hilbert_e=ifft(complex(imag(DFT(1:(length(DFT)/2))),real(DFT(1:(length(DFT)/2)))));
h_e=sqrt(e_downsample(:).^2-hilbert_e(:).^2);
[pks,locs] = findpeaks(abs(h_e));
figure,
plot(abs(h_e));

I will really appreciate, if anyone help me in cases 3 and 4.

As far as I understand, you have a periodic signal, or a quasi periodic signal. And as far as I understand, you have lots of false positives when you detect the peaks.

You need a better way to identify peaks, try this  one...


How can a zero-padded length n signal be truncated to n coefficients from which a signal can be reconstructed? A wavelet transform is defined for infinite length signals. Finite length signals must be extended in some way before they can be transformed. I know that periodic replication and zero padding are appropriate for signals that begin and end on the baseline, while mirror -image replication and linear extrapolation provide continuity at the boundaries for signals that do not begin or end on the baseline. Periodic replication either wraps around or reflect the signal detail in the region beyond the boundaries and this distorts the interpretation of the transform coefficients near the boundaries.

I have a time series of limited duration not extending beyond the signal range and that is not a power of two as required by the transform (using R package wavethresh, function wst(), packet-ordered non-decimated wavelet transform). Zero padding seems to be the only way forward for signals that begin and end on the baseline. Also, zero padding makes no assumptions of the signal after the boundaries describing only the signal, however zero padding results in a non length preserving transform (one in which the transform vector is longer than the signal vector) and large perturbations in the transform space are not reflected in the signal space.

By doing zero padding (added at the beginning and at the end of the series) my question is how can I truncate the zero-padded signal of n coefficients to obtain n coefficients as be able to reconstruct the signal exactly.

From a theory standpoint, that ought to be possible. It's not bound to be simple, however. Zero padding will scale and interpolate the FD representation of your signal. If you remove coefficients the right way, you should be able to reconstruct what's missing.

Let's say you doubled the signal length by padding. You can then remove every other bin of the fourier transformed (starting with the second) and get exactly the result you would have without padding.


What scale should a detected SIFT point be extracted at? I've heard this from a few different researchers but I can't find a definitive value.  Given the "scale" extracted by David Lowe's SIFT detector, we need to extract a descriptor at a specific scale.

Is it 6 times the detector "scale"?  3?  1?  I've heard conflicting values.  What is known about this?

If I understand correctly what you are asking --
In general, the feature is found at the same scale as the SIFT detector says, but
in David Lowe's SIFT, the image is pre-smoothed with sigma: 0.5, so, you need to "subtract" this amount of smoothing from the sigma, so the "real" scale could be:

sqrt(sigma^2 - 0.5^2)

where sigma is the scale the feature was found at.
you have an excellent implementation of SIFT detector and descriptor at: www.vlfeat.org


Generate and Recognize Graphical Patterns I'd like to develop an algorithm for recognizing and generating graphical designs. As an example, consider these examples from one class:

class 1 example 2 http://185.49.85.137/Kave/up/Decorative-Leaves.jpg

Two examples from another class:

I want to detect which class the query image belong to, as well as possibly generate a new sample from the class.

Do you have any idea how this could be done?

I think deep learning comes very handy to solve such problems in these days. Simply create a convolutional neural network (CNN) architecture. Your input are the images, which are resized to one common dimension. The output is the class (type) of the image. For a CNN, this problem is rather simple and if you don't have a lot of training examples, a simple architecture such as LeNet would suffice. Of course, if you have a lot of classes, maybe a larger network (such as AlexNet) might be more appropriate. In either case, you could benefit from a large set of libraries, which are mostly coded in Python. To name a few: Tiny-DNN, Caffe, Tensorflow and Torch.


Design of efficient digital interpolation filter I came across this paper entitled "Design of Efficient Digital Interpolation Filters and Sigma-Delta Modulator for Audio DAC" where the author oversamples an input frequency, fsig = 1kHz with ratio L = 128 and update frequency, fsi = 64kHz. The interpolation filter specification is given by:

passband ripple = 0.001dB for frequency <0.45*fsi and
stopband attenuation = 174dB for frequency > 0.55*fsi.

For two-stage system using IIR Butterworth:

L1 = 2 and L2 = 64, how to get filter order of N1 = 77 and N2 = 21 respectively?
L1 = 16 and L2 = 8, how to get filter order of N1 = 120 and N2 = 7 respectively?

Butterworth low-pass filters are not going to work for this, as demonstrated in the following. You'd need to use other types of filters, for example linear-phase finite impulse response (FIR).

The magnitude frequency response of an order $N$ discrete-time Butterworth low-pass filter can be approximated by that of an analog prototype (from Wikipedia, with a small notation change):

  $$G(\omega) = {\frac{1}{\sqrt{1+{\omega}^{2N}}}},$$

with normalized frequency $\omega$. The Butterworth filter magnitude frequency response is strictly decreasing. For a pass-band edge gain of $a = 10^{-0.001/20}$, that is, -0.001 dB, we'd employ a further gain factor of 0.001 dB to have a gain of 0.001 dB, that is, $\frac{1}{a}$ at 0 Hz. With this gain factor the magnitude frequency response is $\frac{1}{a}G(\omega)$. Then the transition band begins at the frequency $\omega_0$ found by solving:

$$\frac{1}{a}G(\omega_0) = a$$
$$\Rightarrow \omega_0 = \left(\frac{1 - a^4}{a^4}\right)^{\textstyle\frac{1}{2N}}$$

Then we plot the magnitude frequency response at the beginning of the stop band (the end of the transition band) at frequency $\omega_1 = \frac{0.55}{0.45}\omega_0$ as function of filter order $N$:

Figure 1. $\frac{1}{a}G(\omega_1)$ as function of $N$.

Based on the graph, -174 dB of stop band corner gain would require a filter order of about 120. I don't think such a high order is feasible for a recursive filter as it will lead to numerical problems. Filter types other than Butterworth will be able to meet the specification with a lower filter order.


How to filter key clicks? I am the author of an amateur radio application that produces a waterfall display across a range of frequencies. Time is on the $x$-axis, frequency is on the $y$-axis, and the relative strength of each signal is depicted by the intensity of color.

Some of the signals have, what are called, "key clicks". That is, at the start and end of each signal, the bandwidth of the signal is wider than the steady-state bandwidth. An example is shown below:

[1] http://www.kkn.net/~n2ic/key_clicks.jpg

These are actually present in the signal - They are not artifacts of my DSP processing. I would appreciate suggestions on how to filter out the key clicks.

Assuming from your screenshot that the frequency magnitudes of "key clicks" are always higher (brighter) than the normal signal, you could just employ a basic Limiter? In other words, if a frequency magnitude value exceeds a certain threshold, set that frequency value to its last (valid) value. If you wanted to keep the plot from looking unnatural, you could simulate the data where the clicks are or use data from previous samples.


FIR design : compute coefficients for the same frequency response at another sampling frequency I would like to know if there's a standard way to obtain the FIR coefficients for a sampling frequency $F'$ when the coefficients are for a sampling frequency $F$.

I would like to implement the filter described in ITU R BS.1770-3
(the oversampling FIR filter, page 18).
The filter coefficients are given for a sampling frequency of 48 kHz.

What would be the way to derive the coefficients for another sampling frequency?

Edit december 28 2012

I wonder if there's a meaning in my request.
The described filter is used on a 48 kHz sampling frequency audio signal that has been transformed into a 192 kHz sampling frequency signal by stuffing 3 zeros between each original sample. It is the second step of an oversampling method.
The signal $[x_0, x_1, ..., x_n]$ is first transformed into $[x_0, 0,0,0, x_1, 0,0,0, ...,x_n,0,0,0]$. The filter I have mentioned is applied on this zero-stuffed signal.

If I understand correctly (which might not be the case), the cut-off frequency of this low pass filter can be expressed relative to the sampling frequency.
Given this context, and the desired result which is a sort of evaluation of the reconstruction of the audio signal into the analog domain, wouldn't this relative to sampling frequency cut-off frequency be appropriate for any sampling frequency ?

Let's say the cut-off frequency of the filter is $\alpha \times F_s$ ($\alpha$ should be around 0.125 for a 22 kHz analog upper bandwith ?), wouln't this $\alpha$ value be the good one for others sampling frequencies ?
If I have an original signal at 8 kHz sampling frequency, a four time oversample would lead to a 32 kHz sampling frequency, and I can suppose that the reconstructed analog signal would have an upper bandwith limit at approximately 4 kHz ?

End of Edit december 28 2012

There is no direct way of converting filter coefficients between two sample rates while maintaining the exact transfer function (at least where it's properly defined).
Re-sampling the impulse response can be done but will often result in extra latency, a longer filter, and some change in the frequency response.

In this case, however, you can derive the original filter specification from the filter coefficients and than re-design the filter using the same specification at a different sample rate. By visual inspection we can find that the ITU filter is an equiripple filter with

pass band edge: 20 kHz
stop band edge: 28 kHz
pass band ripple: 0.1dB
stop band attenuation: 40 dB

Let's say you want that filter at 4*44.1kHz instead of 4*48kHz. In Matlab you would do the following:

%% match the filter at 44.1 kHz
fs = 4*44100;
d = fdesign.lowpass('Fp,Fst,Ap,Ast',20000,28000,0.1,40,fs);
hd = design(d,'equiripple');
% convert to polyphase
h4 = 4*reshape(hd.Numerator',4,12)';

and the coefficients would be

    0.0057292     0.001552   -0.0015723   -0.0049864
   -0.0056096   -0.0017098    0.0050592    0.0099256
    0.0080942   -0.0011323    -0.012505    -0.017267
   -0.0093558    0.0088501     0.025828     0.027239
    0.0063161    -0.028386    -0.054115    -0.045155
     0.011432      0.10497      0.20287      0.26531
      0.26531      0.20287      0.10497     0.011432
    -0.045155    -0.054115    -0.028386    0.0063161
     0.027239     0.025828    0.0088501   -0.0093558
    -0.017267    -0.012505   -0.0011323    0.0080942
    0.0099256    0.0050592   -0.0017098   -0.0056096
   -0.0049864   -0.0015723     0.001552    0.0057292

However, there are a bunch of potential problems:

In this case we simply lucked out that the number of coefficients for the new sample rate was 48 as well. In other cases you may have to adjust the stop band frequency to dial in the desired number of coefficients.
This particular filter specification has been derived from the original sample rate, i.e. pass band and stop band are symmetrically spaced around the Nyquist frequency of 24 kHz. That is intentional. If your sample rate is different, you need to understand WHY it is different and whether the original specification is still appropriate. In the case of 44.1 kHz you may consider placing pass and stop band symmetrically around 22.050 kHz and maybe also consider lowering the pass band and/or increasing the number of coefficients.


How to Flip Spectrum Around DC? Is it possible to flip a signal's spectrum around DC?  I have a simple spectrum that I made up (MATLAB code):

spectrum =   [-1+4i 0+3i 1+2i 2+1i 3+0i 4-1i 5-2i 6-3i 7-4i 8-5i]
timeDomain = ifft(spectrum);

I looked on this website (http://www.dsprelated.com/showarticle/51.php) and the given proof states this can be done in 3 ways:

Invert the Q channel (14)
Swap the I and Q channels (15)
Invert the I channel (16)

I tried those 3 ways on my signal:

1)  Invert the Q channel

negQtd   = real(timeDomain) - 1j * imag(timeDomain);
negQSpec = fft(negQtd)

Output:

[-1-4i 8+5i 7+4i 6+3i 5+2i 4+1i 3+0i 2-1i 1-2i 0-3i]

2)  Swap the I and Q channels

swapTD   = imag(timeDomain) + 1j * real(timeDomain);
swapSpec = fft(swapTD)

Output:

[4-1i -5+8i -4+7i -3+6i -2+5i -1+4i -0+3i 1+2i 2+1i 3+0i]

3)  Invert the I channel

negItd   = -1 * real(timeDomain) + 1j * imag(timeDomain);
negISpec = fft(negItd)

Output:

[1+4i -8-5i -7-4i -6-3i -5-2i -4-1i -3-0i -2+1i -1+2i -0+3i]

Each of these does flip the spectrum of the signal, but it also modifies the spectrum:

1) Inverting Time Domain Q channel also negates/inverts the Frequency Domain Imaginary part.

2) Swapping the I and Q channels also swaps the Frequency Domain Real and Imaginary parts.

3) Inverting the Time Domain I channel also negates/inverts the Frequency domain Real part.

Is there some other way to simply flip the spectrum around DC or should I try a two step process?

The article Handling Spectral Inversion in Baseband Processing that you refer to in your question is not about simple inversion of the frequency axis, but it is about inversion of the frequency axis and conjugation in the frequency domain. This is necessary because due to different mixing conventions in several conversion stages, the upper and lower side-bands can be exchanged, and the two side-bands exhibit conjugate symmetry, because the transmitted signal is real-valued. So one side-band can be obtained from the other by inverting the frequency axis and by additional complex conjugation (in the frequency domain).

This process of inversion and conjugation in the frequency domain corresponds to complex conjugation of the complex baseband signal in the time domain, which is equivalent to simply inverting the $Q$ component. Since a phase shift of $180$ degrees is irrelevant, one can equivalently invert the $I$ component. A phase shift of $90$ degrees corresponds to a multiplication with $j$ in the time domain, and if this phase shift can be tolerated, you get the additional option of swapping the $I$ and $Q$ channels.

If you want pure inversion of the frequency axis without conjugation, then you indeed need to invert the time axis of the corresponding time domain signal (as pointed out in a comment by @curiousStudent). Luckily this is usually not needed in practice because this can't be done easily in a real-time system.

As a final note, all three options in your question do flip the magnitudes of the spectra around DC, it is just the phases that behave differently in all three cases.


How to express one image in terms of another one I have two (black and white) images of identical size - let's say 128x128 pixels. I'm interested in expressing Im2 in terms of components of Im1 such that as opposed to specifying Im2, I can specify some values (less than number of pixels in Im2, otherwise what's the point) which, together with the knowledge of Im1 would allow me to construct an approximation of Im2. Let me demonstrate by giving you my first idea (which, spoiler alert, didn't work). In pythonish pseudocode:

# Singular Value Decomposition of Im1 and Im2
U1, S1, V1 = svd(Im1)
# Function to minimise
def fit_function(S):
    Im_fit = U1*diag(S)*V1
    diff = Im2-Im_fit
    return sqrt(mean(diff^2)) # root-mean-square difference
# Find S which fits Im_fit to Im2
Sfit = myLittleMagicSolver(fit_fun)

In other words; use singular value decomposition on Im1 and Im2. Then find numbers other than singular values of Im1 which together with U1 and V1 matrices approximate Im2. This way I could specify only 128 values - the main diagonal of the fitted S matrix.

This, however, does not work. My fitted images are nowhere close to their target of Im2. My question is then, is that a problem someone has tackled already? I am not fussed about the SVD/PCA approach. I just want to be able to express Im2 in terms of Im1.

My intended use is to train a neural network to recognize patterns in my particular dataset where one image is an input to the network, and multiple images are the output. This, of course, makes for a pretty hefty neural net. I, therefore, want to teach the network to recognize the components of the output images in terms of the input instead.

I found an answer which is good-enough for me. As @Stanley Pawlukiewicz has pointed out in the comments, this is hard to do for a general case when there is little correlation between the images. I, however, want to work with real images of actual things in the real world. This means there will be a lot of low frequency components (that's why jpeg compression works). If I now formulate the problem as:
$$
I = U_{s} S  V_{s}'
$$
where $I$ is the matrix representing the image, $U_{s}$ and $V_{s}$ are the matrices of the svd decomposition of the source image, and we are trying to find $S$. Then, the inverse problem (which is the one we actually want to solve) is:
$$
S = U_{s}^{-1}I(V_{s}')^{-1}
$$
This is the easy part (using Matlab, python or otherwise).

I then took a good look at how are the solutions ($S$) to this problem are distributed, and they tend to cluster around the left and upper edge of the S matrix. Therefore, I can approximate S by the subset of S centered around these values. I (slightly arbitrarily) have chosen to use a quarter of a superellipse to set the cutoff value with a nice constant to play with. I can get S down to 10% of its elements and retain most of the broad features of the image, which is all I wanted.

I've put my code up online so you can check it out. Hope this helps someone. :)


Fourier Transform gives unexpected results: signal reversal and conjugation As far as I understand the math you can reverse a real valued signal by Fourier transformation, taking the complex conjugate of the result and inverse Fourier transformation, i.e.
\begin{equation}
\tag{1}
f(-x) = \mathcal{F}^{-1}(\overline{\mathcal{F}(f(x))})
\end{equation}
where $\mathcal{F}$ is the fourier transformation operation and the overline denotes complex conjugation. This all makes sense and from a mathematical perspective this isn't very useful. However, from a DSP perspective its very useful because the autocorrelation of a signal can be conveniently cast as the convolution of a signal with a time-reversed version of itself. If you're using an FFT convolution it means that you only have to calculate one FFT of the signal.

When I tried to confirm that equation (1) was working in python I found that the result is shifted by one sample. Here's the example code:

import numpy as np
from numpy.fft import rfft, irfft
x = np.arange(128)
y = np.sin(2 * np.pi * x / 128 * 4 + np.pi / 3)
y_rev = y[::-1]
y_rev_fft = irfft(conj(rfft(y)))
plt.plot(y, label="Original Signal")
plt.plot(y_rev, label="Straight Reverse")
plt.plot(y_rev_fft, label="FFT Reverse")

Is this the expected result? If so, can the conjugate still be used in the FFT correlation algorithm I described?

Also, is this question better suited to another stackexchange site?

Try plotting the signal and the straight reversed signal in points rather than lines and you will see that, unlike with reversal via FFT, none of the values stay the same throughout the reversal. The straight reversal flips the 128-long data around index 63.5 or equivalently around index -0.5 if you think of the data as periodic with a period of 128, rather than around index 0 (or 64) as in reversal via FFT.

According to user David Hoffman's tests both

roll(irfft(conj(rfft(y))*rfft(y)), -1)

and

irfft(rfft(y[::-1])*rfft(y))

work for autocorrelation, roll being numpy circular shift.


Why ignores constant when applying integral into differential-equation.? I studied how to draw 'Direct-form I' and 'Direct-form II' from differential equation.

The book taught me like this:

$$\mbox{Give equation: }{a_0}{y(t)}+{a_1}{y'(t)}+{a_2}{y''(t)} = {b_0}{x(t)}+{b_1}{x'(t)}+{b_2}{x''(t)}$$

step 1. change from differentiator into integrator, because integrator is cheaper than differentiator.

\begin{array}{}
{a_0}{y^{(2)}(t)}+{a_1}{y^{(1)}(t)}+{a_2}{y(t)} = {b_0}{x^{(2)}(t)}+{b_1}{x^{(1)}(t)}+{b_2}{x(t)} &
\begin{cases}
x^{(i)} \mbox{ is i-times integral of x}\\
y^{(i)} \mbox{ is i-times integral of y}
\end{cases}
\end{array}

step 2. change form like this: $y(t)=...$

\begin{array}{}
\displaystyle y(t)=\frac{1}{a_2}\left( {b_0}{x^{(2)}(t)}+{b_1}{x^{(1)}(t)}+{b_2}{x(t)}-{a_0}{y^{(2)}(t)}-{a_1}{y^{(1)}(t)}\right)
\end{array}

step 3. draw

My question is why doesn't consider constant value in step 1?

The constants would be the initial conditions, such as an initially charged capacitor or an inductor with an initial current. If all constants are set to zero it means that you have zero initial conditions, which is usually the standard assumption unless stated otherwise.


Multiple parallel IIR bandpass filters with different centre frequencies combined into one filter for optimisation... possible? Before I begin, I have read these already:

IIR filter parallelization
Concept of combining multiple FIR Filters into 1 FIR filter

I am processing audio using Csound in a high performance real time environment. One part of my code uses 10 butterworth bandpass filters in parallel, each with a different centre frequency and bandwidth. However, the centre frequencies and bandwidths are all relative to each other. The centre frequency of each filter are set whole number multiples or divisions of the main filter's centre frequency (ie, if main filter has cf of 100Hz, the other filters cfs would be 200, 300, 400, 500, 600, 50, 33.333, 25, and 20). The bandwidths of each filter are all a specific division of the associated centre frequency (ie, cf 100 = bw 10, cf 200 = bw 20, etc).

The main centre frequency is not fixed, but seeing as the filter centre frequencies and bandwidths are all relative to the main filters centre frequency and bandwidth, could I design a filter that combines all these filters into one? And most importantly, if I could, would it make any major difference to performance?

My best guess as to how to achieve this would be to combine the filter equations into one. I don't know if this would actually work, but if it did, I don't feel like it would improve the performance by much.

EDIT:

The sample being filtered is rich in harmonic and in-harmonic frequencies. It is borderline noisy, so the sum of the outputs of the BP filters is quite interesting, and the harmonies aren't always perfectly related. This quality I would like to keep.

So, what I took away from my multirate lectures is:

You don't use IIRs in that case, because the recursive part isn't shareable among multiple filter implementations (being specific to the output of each filter), and hence, you end up with a lot more multipliers than you'd have if you went for a FIR system and removed redundancies. This is especially true for hardware/FPGA implementation, where placing a wire to use a result twice is cheap. You're doing software, so saving a few registers worth of filter taps really won't pay -- you're breaking the linear nature of memory access that way, and that will impose a performance hit that is much, much worse.

In fact, I'd challenge you to look at how you can understand your filter bank as an array of filterbanks with identical filters.

For example, assume that all your BPFs have equally spaced center frequencies $\frac{f_\text{sample}}{N}$. Also, they might have different bandwidths, but since you say they are related to each other, let's assume there's one "maximum bandwidth" that is a multiple of all smaller bandwidths.

Now, implementing a polyphase filterbank ("channelizer") with a BPF of that maximum bandwidth and the constant frequency spacing is very efficient. Basically, it's like taking a single very BPF at full rate, add one FFT, and you get the complexity of the system for up to $N$ channels (note: this is $N-1$ free channels!). Then, take the channels whose bandwidth you'd need to be reduced further. Now, assume these $m<N$ channels all need to be reduced in bandwidth by a factor of $4$, for example. So that's $m$ identical quarterband lowpass filters! Implementing identical filters in hardware is pretty efficient, usually, especially because after the first filterbank, sampling rate might already be reduced (so you can multiplex multipliers); now, you're seemingly doing audio on a PC:

Take these $m$ channels, and just implement individual filters. There's nothing much you can or need to optimize here; assuming this is on something that has a PC-like CPU, caching will keep your filter taps close to CPU registers, and and only the actual samples will need to go through memory buses. Also, not that things like $\frac{1}{M}$ low pass filters are relatively relaxed in transition width requirements, and hence, FIR implementations are short in taps, and thus fast.

As usual in signal processing, decimate whenever possible! The polyphase filterbank might already decimate, but also your $m$ channel low pass filters should do that, too, and in a clever manner (a polyphase lowpass decimator runs at the output, not the input sampling rate!). Not processing samples you don't need is probably reason Nr. 1 that modern DSP is so capable.


Variational Regularization Method in Image Processing I would like to understand better variational regularisation methods in image processing. In particular, the formulas in this image:

Why formula (3.13)? In the notes that I read I cannot find anything about its background, why is the solution of (3.12) given by (3.13). Why do we need this part $$\alpha||f||^2$$?

This is an example of the Fidelity Term and Prior Term model.
In many Inverse Problems we assume some model on the additive noise. This part is modeled by the Fidelity Term ($ \mathcal{D} \left(A \boldsymbol{f}, \tilde{\boldsymbol{g}} \right) $ in your example). For Gaussian Noise it is given by Least Squares Term:
$$ \frac{1}{2} {\left\| A \boldsymbol{f} - \tilde{\boldsymbol{g}} \right\|}_{2}^{2} $$
The prior terms is usually something else we know about the data in question. It can be something about its distribution, about its derivatives or in the case of Variational Models something about its response to some operator (For example, Derivative).
In 3.12 above they use the Tikhonov Model which I wrote about in your other question.


Real signal from I,Q, sample rate and center frequency I am given a recording of an ionosonde sounding. The signal $S$ is represented as a series of  I/Q pairs at a sample rate $F_{sr} = 10\,\text{MHz}$.  I am told that the center frequency of the recording is $F_{center} = 7\,\text{MHz}$.

I have a $T=200\,\text s$ recording as $I=[I[0],..,I[T\times F_{sr}-1]]$, $Q=[Q[0],..,Q[T\times F_{sr}-1]]$.

The actual signal $R[t]$ being recorded can vary from $2\,\text{MHz}$ to $12\,\text{MHz}$, that is from ${F_{sr}\over 2} -F_{center}$ to ${F_{sr}\over 2}+F_{center}$.

Let $S=I+Q j$ be the complex representation of the recording.

How do I recover the original signal $R[t]$ at a sample rate of $F_{up}=24\,\text{MHz}$, which is high enough to represent a $12\,\text{MHz}$ signal directly, as a function  $R(S,F_{center},F_{sr})$?

That is, $F_{up}=2F_{sr}+ 2  F_{center}$, so $R=R[0],...,R[T\times F_{up}-1]$.

I gather you are trying to reconstruction a real (non-imaginary) version of your I & Q signal. The following scheme should work:

Up sample your IQ from 10 MHz to 24 MHz. This involves zero padding and applying a low-pass filter to remove the extra images that appear due to the zero padding. See below for more details.
Mix the signal with $\exp(j2\pi7\times 10^6t)$. This shifts the signal to the 7 MHz.
Take the Real component of the result i.e. throw away the imaginary component.

For the upsampling process you have a couple of options. You can upsample by 12 and then decimate by 5. You would up sample by 12, apply the low-pass filter, and then down sample by 5. The filter would need a pass band from 0 - 4MHz (exact upper frequency depends on the nature of your IQ signal). The stop band would need to be at 5 MHz. You might be able to get away with a higher stopband frequency - e.g. if you have no frequency content between 4 MHz and 5 MHz, you could use a stopband starting at 6 MHz. The filter design is slightly dependent on the frequency analysis of your IQ signal and you would have to figure out your requirements for the filter specifications e.g. pass-band ripple and stopband attenuation.

For practical purposes you probably want to break the up sample by 12 into two separate steps e.g. Up sample by 4 and then up sample by 3 - each step would have its own filter applied.

An alternative method is to up sample by a factor of 3 or 4, apply the low-pass filter and then use a spline filter to interpolate the signal to required 24MHz sampling rate.


Output of KISS FFT interpretation I'm receiving PCM data trough socket connection in packets containing 320 samples. Sample rate of sound is 8000 samples per second. I am doing with it something like this:

int size = 160 * 2;//160;
int isinverse = 1;
kiss_fft_scalar zero;
memset(&zero,0,sizeof(zero));
kiss_fft_cpx fft_in[size];
kiss_fft_cpx fft_out[size];
kiss_fft_cpx fft_reconstructed[size];

kiss_fftr_cfg fft = kiss_fftr_alloc(size*2 ,0 ,0,0);
kiss_fftr_cfg ifft = kiss_fftr_alloc(size*2,isinverse,0,0);

for (int i = 0; i < size; i++) {
    fft_in[i].r = zero;
    fft_in[i].i = zero;
    fft_out[i].r = zero;
    fft_out[i].i = zero;
    fft_reconstructed[i].r = zero;
    fft_reconstructed[i].i = zero;
}

// got my data through socket connection

for (int i = 0; i < size; i++) {
     // samples are type of short
     fft_in[i].r = samples[i];
     fft_in[i].i = zero;
     fft_out[i].r = zero;
     fft_out[i].i = zero;
 }

 kiss_fftr(fft, (kiss_fft_scalar*) fft_in, fft_out);
 kiss_fftri(ifft, fft_out, (kiss_fft_scalar*)fft_reconstructed);

 // lets normalize samples
 for (int i = 0; i < size; i++) {
     samples[i] = rint(fft_reconstructed[i].r/(size*2));
 }

After that I fill OpenAL buffers and play them. Everything works just fine but I would like to do some filtering of audio between kiss_fftr and kiss_fftri. Starting point as I think for this is to convert sound from time domain to frequency domain, but I don't really understand what kind of data I'm receiving from kiss_fftr function. What information is stored in each of those complex number, what its real and imaginary part can tell me about frequency. And I don't know which frequencies are covered (what frequency span) in fft_out - which indexes corresponds to which frequencies.

I am total newbie in signal processing and Fourier transform topics.

Any help?

What you might want to investigate is FFT fast convolution using overlap-add or overlap-save methods. You will need to expand the length of each FFT by the length of the impulse of your desired filter. This is because (1) FFT/IFFT convolution is circular, and (2) each index in the FFT array result corresponds to almost all frequencies (a Sinc shaped response), not just one (even if mostly concentrated near one), so any single bin modification will leak throughout the entire frequency response (except certain exact periodic frequencies).


Confusion regarding calculation of fundamental period? I am reading signal processing first and in chapter 3 ex3.8 i came across an example of fundamental period as shown  in attached photo

It apparently shows that signal $$x(t)=\cos^2(4\pi t)$$ has period 0.5 but then it also writes that fundamental period is 0.25

How is he doing that?

Also will be fundamental period if $$x(t)=\cos^n(4\pi t)$$ where n can be 3 or 4 or 5

This seems like more of a semantics problem.

A signal is periodic with time $T$ if

$$x(t+n\cdot T) = x(t), n \in \mathbb{Z}$$

So the signal is periodic in $0.5$ since the for $T = 0.5 \cdot n$ the argument of the cosine is an integer multiple of $2 \pi$. Since it's periodic in $0.5$ it's also periodic in all integer multiples of $0.5$, i.e $1$, $1.5$, $2$ etc.

In this case it's also periodic in $0.25$ since
$$ \cos^2(4 \cdot \pi \cdot t ) = 0.5 \cdot (1+\cos(8 \cdot \pi \cdot t))$$

So any periodic signal has an infinite number of periods, the fundamental one is the smallest one and all the others are integer multiples of the fundamental.


Preserving integral through downsampling I've got a high-resolution image segmentation, which I'm downsampling (and non-linearly deforming) to fit on an ROI image. I typically blur the high-res image before downsampling with a Gaussian-fwhm related to the downsampling factor. The result is high-res and low-res images with the same mean.

I'd like now for my high-res and low-res images to have the same sum. If I scale up the low-res image by the downsampling factor, will local sums be preserved across the image? Is there a better way to downsample and preserving local sums?

Increase the gain of your filter by the downsampling factor.

You do this by multiplying the filter taps by the downsampling factor.  This is equivalent to doing exactly what you have done before and then multiplying the results by the downsampling factor.


Effect of Resizing (Up Scaling) an Image on the Ability to Detect / Recognize Objects I work in a research lab and am currently setting up a system to image objects approximately 30 microns in diameter with a 250x USB microscope, then use image subtraction to isolate the objects as they move past the microscope FOV, and use blob detection to detect the object's presence in the binary subtracted image.

The microscope can capture still images at 1280x720, but video can only go up to 640x480. Ideally I wanted to take rapid still images (~10 FPS) at 1280x720, but various issues are preventing me from doing that. Capturing rapid images at 640x480 on the other hand is easily doable, but I'm not sure if the images would be clear enough for the image processing to resolve the objects at that resolution.

Would capturing the images at 640x480, then resizing them to 1280x720 before processing them, improve the ability to detect/resolve the objects in the image?

It is theoretically impossible to add information by increasing the resolution of an existing image.

So from that point of view -- no.

However, I certainly find it conceivable that if you're using some canned algorithms they may be sensitive to pixelation.  Unless they're really lame, I would only expect this to be an issue if your original images are very sharp -- if the un-differenced objects tend to have blurry edges, either because of the limitations of your microscope or the properties of the objects themselves, then I'd think about trying what you suggest.

Otherwise -- you can't get something for nothing!


Expected value of coupled signals is zero Suppose there are two vector signals, $u$ and $v$, which are functions of time. Given, $E[uv^T]=0$, what is the physical significance of this statement. For clarity, assume $u$ and $v$ are noises in a system like in kalman filter, control noise and measurement noise. E is the expected value and assume that the vector dimensions are entirely compatible.

It means u and v are orthogonal.

If u or v or both have zero mean then they are also uncorrelated.

If they are uncorrelated and jointly gaussian then they are also independent.


Image Noise Detection in MATLAB I am getting error while implementing a simple idea for noise detection in an image. Let  I  denote  noisy  image  and  for  each  pixel of  I  is represented by I(x,y). A Sliding window of size 3X3 centered at I(x, y) is defined. At first we estimate that the center pixel of window is noise or not. To differ the noise and signal, we calculate the mean (m) and standard deviation (s) of the window. If the pixel is between m-s and m+s then it is a signal otherwise noise and we will apply median Filter.

Algorithm:

1.  Read an Image I.
2.  Take a Sliding Window or Mask of size 3X3.
3.  Calculate the Mean (m) and Standard deviation (s) from the Mask.
4.  Calculate the threshold as follows
         t1 = m-s and t2 = m+s

5.  IF     t1 <= I(x,y) and  I(x,y)  <= t2 THEN

           Result(x,y) = I(x,y)
    ELSE
        Result(x, y) = medfilt2(I(x,y)

6.  Repeat the step 3 , 4 and 5 on Entire Image.
7.  Display the Resultant Image.

I tried to implement this algorithm as follows

clc;
close all;
clear all;

I=imread('lena.jpg');
I=rgb2gray(I);
J=imresize(I, [128 128]);
L = imnoise(J,'speckle',0.01);

[m,n]=size(L);
% RETAINING THE CORNER PIXELS/ we can do zero padding also, suggest me which one is good

for i= 1:m+4
    for j = 1:n+4
        L(:,1:2)=J(:,1:2);         % retain the first two columns
        L(:,(m-1):m)=J(:,(m-1):m); % retain the last two columns
        L(1:2,:)=J(1:2,:);         % retain the first two rows
        L((n-1):n,:)=J((n-1):n,:); % retain the last two rows
    end
end

for i=1:m-2
    for j=1:n-2

      P=L(i:(i+2),j:(j+2));
      s=std2(P(:)); % Calculating Standard deviation
      m=mean2(P(:)); % Calulating Mean
      t1=m-s; % Threshold 1
      t2=m+s; % Threshold 2

      if (P(2,2)<=t1 && P(2,2)>=t2)
       iout(i,j)=medfilt2(P(:))
      else
       iout(i,j)=P(2,2)
      end

      end
end
L
iout
imshow(iout, []);
imshowpair(L,iout,'montage');

It runs but there is no change between input and output images.. Image is not filtered.
Please help me.

If you would like to do this faster, what you could try is the neighborhood operations. Take for example the sample code below:

I = imread('tire.tif');
f = @(x) max(x(:));
I2 = nlfilter(I,[3 3],f);
imshow(I);
figure, imshow(I2);

This applies a maximum filtering. You could consider changing your neighborhood function f to your custom algorithm (inside your loop) and this should give you a better speed.

Of course there are a lot of ways of getting speed. Use of integral images, moving averages and precomputations might save you from the burden of heavy block processing. But this would be way a lot to write in here. For more details check out these papers:

http://www.csse.uwa.edu.au/~shafait/papers/Shafait-efficient-binarization-SPIE08.pdf
http://arxiv.org/ftp/arxiv/papers/1201/1201.5227.pdf

Cheers,


What's the formula for choosing FFT sample rate and window length to make sure the frequency I want is in the center of a bin? I don't know why I am having such a hard time with this math, but I am. I am doing phase-sensitive detection of a few sensors, which are triggered by a modulated light source. The modulation frequency f0 changes by up to several Hz over the course of the day based on ambient conditions, and is typically around 1370-1390 Hz. I'm doing an FFT to extract the sensor amplitudes at f0. With my system, I'm limited to an upper sampling rate of 10 kHz for up to 2 seconds, but no longer.

Given that I can dynamically change f0, sample rate, and window length, what is the formula for making sure that f0 is in the center of a bin?

The concept you need is to have a whole number of cycles within your sample frame.  The cycle count will be the bin index.

For the math, consider carrying units:

$$ \frac{Cycles}{Frame} = \frac{Cycles}{Second} \cdot \frac{Seconds}{Frame} $$

Where:

$$ \frac{Seconds}{Frame} = \frac{  \frac{Samples}{Frame} }{ \frac{Samples}{Second} } $$

Solving for your sampling Rate:

$$ \frac{Samples}{Second} = \frac{  \frac{Samples}{Frame} }{ \frac{Seconds}{Frame} } $$

$$ \frac{Samples}{Second} =  \frac{ \frac{Samples}{Frame} \cdot \frac{Cycles}{Second}  }{ \frac{Cycles}{Frame} }  $$

If you have a pure tone, there isn't a need to center it on a bin to find your parameters, although it is much easier to do it that way.

To find the frequency, you can use the formulas in these blog articles:

Exact Frequency Formula for a Pure Real Tone in a DFT
Two Bin Exact Frequency Formulas for a Pure Real Tone in a DFT
Improved Three Bin Exact Frequency Formula for a Pure Real Tone in a DFT

Once you have the frequency, you can calculate the amplitude and phase with the formulas in this article:

Phase and Amplitude Calculation for a Pure Real Tone in a DFT: Method 1

You should also realize that you don't need to call a full FFT.  You simply need to calculate the two or three bin values you are interested in which is generally a lot quicker.


Amplitude Estimating Using a Windowed DFT Let's say we want to estimate the amplitude A of a mono-frequent signal using a windowed DFT. The frequency of the signal is unknown, and the frequency resolution of the DFT is limited, thus it cannot be expected that there is a DFT bin that matches this frequency directly. Therefore, A is estimated from the amplitude of the next closest DFT bin. Now the question is,

should the frequency response of the window exhibit a narrow or a wide mainlobe to obtain a good estimate, and why?

I have several blog articles that solve this problem exactly in a theoretical sense, and quite accurately in an implementattion.

3. DFT Pure Tone Frequency Formulas

Exact Frequency Formula for a Pure Real Tone in a DFT
Two Bin Exact Frequency Formulas for a Pure Real Tone in a DFT
Improved Three Bin Exact Frequency Formula for a Pure Real Tone in a DFT
A Two Bin Exact Frequency Formula for a Pure Complex Tone in a DFT
Three Bin Exact Frequency Formulas for a Pure Complex Tone in a DFT

4. DFT Phase and Amplitude

Phase and Amplitude Calculation for a Pure Real Tone in a DFT: Method 1
Phase and Amplitude Calculation for a Pure Complex Tone in a DFT
Phase and Amplitude Calculation for a Pure Complex Tone in a DFT using Multiple Bins

The rest of the articles are pretty good too. ;-)

I posted one this morning:

Angle Addition Formulas from Euler's Formula


Is it possible to reduce the complexity of radix-2 FFT if the input vector contains identical elements? my question is about reducing the complexity of radix-2 FFT when the input vector has a specific structure.

For an input vector of x with N elements, the complexity is given by O(N log2 N).

My input vector y with N elements has the structure of:
y = [x1 x2]T

where,

x1 = x2

and x1 and x2 vectors have N/2 elements. Is it possible to reduce the complexity of radix-2 FFT for this y vector? I think it is possible to not use one of the N/2 DFT blocks of the FFT algorithm in this case. Can anyone show me what is the complexity for this vector? Thanks in advance.

if

$$\text{FFT}\{[a, b, c, d]\} = [p, q, r, s]$$

then

$$\text{FFT}\{[a, b, c, d, a, b, c, d]\} = [p, 0, q, 0, r, 0, s, 0]$$

or
$$=[2p, 0, 2q, 0, 2r, 0, 2s, 0]\text,$$

depending on normalization.

Edit:
Using this property, we can make FFT nearly 2x faster on this specialized data. But the complexity is still $O(n \log n)$. I mean, big O notation doesn't count constant factor.


M-point FFT Amplitude of N-samples of a sinusoid Suppose you have N samples of a sinusoid of frequency $f_0$.

$$
x[n] = A\sin[2\pi (f_0/f_s)n ]
$$

Suppose you take an M-point FFT (where M may be greater than N) of your N samples of x. What will the amplitude of the peak of the spectrum be?

Reference: http://www.designnews.com/author.asp?section_id=1419&doc_id=236273&print=yes

In general the FFT of a sinusoid function is NOT a delta function. The FFT's maximum will be around the sine frequency but the exact value will be dependent on how well the frequency of the sine lines up with the FFT grid and there will be non-zero values for most other frequencies.

You only get a delta function if the sine frequency is an integer multiple of your FFT bin spacing (= sample rate divided by FFT length). In this case the amplitude will N/2 (for conventional FFT scaling) for both the positive and the negative bin.

If your FFT length is shorter than the sine wave you have to truncate and the above considerations apply. If the FFT length is longer, you will probably zero pad. This will often result in substantial increase of the energy at the "other" frequencies and the effects depend heavily on the details of the zero padding.

In all cases Parseval's Theorem holds: total energy in the time domain equals total energy in the frequency domain. That's how you can derive an amplitude of N/2 for the special case that the sine wave frequency falls on an FFT bin: A sine wave of amplitude 1 has an RMS value of sqrt(.5). A buffer of N samples has a total energy of sqrt(.5)*N. This is the total energy. In the frequency domain this gets split between the bins at +f and at -f. Each bin receives have of the total energy which is sqrt(.5)*sqrt(.5)*N = N/2.


Analytical design of a bandstop Butterworth filter I am trying to analytically design a bandstop Butterworth filter. I am trying to use:
$$|H(s)|^2 = \frac{1}{1+\big(\frac{s}{j\omega_1}\big)^{2N}}+\frac{1}{1+\big(\frac{s}{j\omega_2}\big)^{-2N}}$$
to derive a transfer function. I would like to know if I'm correct and how to proceed from here.

it looks like you are designing a parallel highpass and lowpass filter.  i don't think that's the best way to do it.  instead, take your butterworth LPF prototype

$$\Big|H(j\Omega)\Big|^2 = \frac{1}{1+\Omega^{2N}}$$

and make this substitution for normalized $\Omega$:

$$ j\Omega = s \leftarrow \frac{Q}{\frac{s}{\Omega_0} + \frac{\Omega_0}{s}} $$

$Q$ is inversely proportional to the bandwidth of the bandstop filter  (if $N=1$, then $Q$ is the actual "Q" of the second-order filter.  $\Omega_0$ is the angular frequency of where the bandstop is tuned to.


Is Smith's book still usable and up-to-date? I'm learning DSP with Smith's book, but I discovered the second edition is from 1999, is this book still usable today?

Do you think that there are a lot of new things have been introduced in the last 10 years? :-)

I like this book. Book is quite accurate and quite pragmatic.  I think this book is a great introduction to DSP. But feel free  to use another books, especially if you want to known about the modern hardware equipment.


Hi, If I really want to focus on the mathematics of signal procesing, Should I do a Masters in EE or Math? I'm an undergraduate electrical engineering student who has an interest in mathematics of signal processing. I want to be able to go to a graduate school where I can focus on more the mathematics and theory behind signal processing. Can this be achieved in a Masters in EE or is it more efficient to go to my Masters in APPLIED MATH?

As stated in the comments, this is more about the content than the label. The best is to play with both fields if you want to invest in the future: signal processing is borrowing  increasing from mathematical tools, for instance convex and non-convex optimization, learning, topological and geometrical tools (for dimension reduction), graphs, etc., that might  not be taught in traditional EE. Such techniques can help with the growing quantity of data with little models and knowledge.

As a personal anecdote, I attended a lecture on "math for EE engineers" for the beginning of a signal processing PhD program in a US university. It was about Euclidean vector spaces (in finite dimension of course). A young PhD student asked whether the basis vectors did belong to the space. I was quite surprised by the level of the question. But not as much as when I hear the teacher say, after a long pause: "mmmm, I'd say... yes, let me think about it" (such a question is indeed interesting in infinite dimension spaces). There are similar  stories from colleagues attending other MsC programs.

So after picking the right topics, I strongly suggest  that you either complement an EE program with maths, to get better insight and intuition about the theoretical    tools, or dive in maths, while practicing a lot on "real DSP problems" aside. The gap is large between mastering  Fourier theory and being able to properly estimate and interpret  of a power spectrum of a noisy signal.
There are plenty of online resources to assist you in that task.

I really appreciate the concept  behind the law of excessive learning in mathematics, stated by Alexander Borovik:

  To be able to use maths at certain level it is necessary to learn it
  at the next level


Detecting corners using structure tensor matrix How does eigenvalues of structure tensor matrix help in detecting corners?
For eg, in Harris corner detection, we consider the matrix
$$M = \sum_{w} w \cdot
\begin{bmatrix}
I_x^2 & I_x I_y \\
I_x I_y & I_y^2
\end{bmatrix}
$$
We say that if both eigenvalues are high, then it is a corner. If one of them is high, then it is an edge. If both are low, then it is neither corner nor edge. How do we come to this conclusion?

I'll simply present a series of results relevant to this question.

Note that eigenvalues of this matrix are non-negative

Let $x = [u\ v]^T$
$$x^TMx = [u\ v]
\begin{bmatrix}
I_x^2 & I_x I_y \\
I_x I_y & I_y^2
\end{bmatrix}
\begin{bmatrix}
u \\
v
\end{bmatrix}
 = u^2 I_x^2 + 2uv I_x I_y + v^2 I_y^2 = (uI_x + vI_y)^2 \geq 0$$
Hence $M$ is positive-semidefinite and hence eigenvalues are always non-negative.

Eigen vector of the matrix corresponding to largest eigen-value gives the direction of maximum change

Consider a symetric matrix $A$. We want the direction of maximum change
i.e.
$$\max_{x} \frac{\Vert Ax \Vert}{\Vert x \Vert}$$
Maximizing the above is equivalent to maximizing $\Vert Ax \Vert$ by fixing $\Vert x \Vert = 1$. Also maximizing $\Vert Ax \Vert$ is equivalent to maximizing $\Vert Ax \Vert^2$
$$ \max_x  \Vert Ax \Vert^2 \quad s.t. \quad \Vert x \Vert = 1 $$
Lagrangian gives
$$ \max_x \quad x^TA^TAx + \lambda (1 - x^Tx) $$
Differentiate w.r.t. $x$ and equate to $0$
$$ 2A^TAx - \lambda x = 0 $$
$$ A^TAx = \lambda x $$
Hence
$$ \Vert Ax \Vert^2 = \lambda \Vert x \Vert^2 = \lambda$$
To maximize this quantity, choose $\lambda$ as the largest eigenvalue and hence $x$ is the eigenvector corresponding to this eigenvalue.
Since $A$ is symmetric $A^TA = A^2$ and hence eigenvalues of $A$ are sqrt of eigenvalues of $A^TA$ and eigenvectors remain same.

Eigenvectors of $M$ are orthogonal
Since $M$ is symmetric, it's eigenvectors are orthogonal. See here.
We found that the eigenvector that corresponds to largest eigenvalue is in the direction of maximum change. Hence the second eigenvector will be in the perpendicular direction to this eigenvector.

Eigenvalues in smooth region
In smooth region, $I_x=I_y=0$ and hence eigenvalues are $0$

Eigenvalues at edges
Let's say the edge is oriented at an angle of $\theta$.
Hence
$$I_x=c \ \sin(\theta);\quad I_y=c \ \cos(\theta)$$
$$ M = c^2
\begin{bmatrix}
\sin^2(\theta) & \sin(\theta)\ \cos(\theta) \\
\sin(\theta)\ \cos(\theta) & \cos^2(\theta)
\end{bmatrix}
 $$
Eigenvalues are given by
$$ (\sin^2(\theta) - \lambda)(\cos^2(\theta) - \lambda) - \sin^2(\theta)\cos^2(\theta) = 0 $$
$$ \lambda^2 - 2\lambda + \sin^2(\theta)\cos^2(\theta) - \sin^2(\theta)\cos^2(\theta) = 0 $$
$$ \lambda^2 - 2\lambda = 0 $$
$$ \lambda = 0, 2 $$
Hence at an edge, one eigenvalue is $0$ and the other isn't.

Eigenvectors at edges
Consider a vertical edge. Hence
$$I_x=c;\quad I_y=0 $$
$$ M =
\begin{bmatrix}
c^2 & 0 \\
0 & 0
\end{bmatrix}
 $$
Eigen vector corresponding to largest eigenvalue is $v = [1\ 0]^T$ i.e. gives the normal to the edge.

Also note that, for the eigenvector in the direction where there is change in pixel intensity, eigenvalue is high. For the eigenvector in the direction where there is no change in intensity, eigenvalue is low. So, eigenvalues represent the amount of change in pixel intensity observed along the direction of it's corresponding eigenvector.

Eigenvalues and eigenvectors at a corner
Let's consider a corner where a vertical edge crosses a horizontal edge. The first eigen vector gives normal to the stronger edge. The second eigenvector, which is perpendicular to the first gives normal to the other edge. Since the change is high along both directions, the eigenvalues corresponding to both the eigenvectors will be large.

Thus, to summarize

In smooth regions, eigenvalues are close to 0.
At edges, one eigenvalue is large and the other is close to 0.
At corners, both eigenvalues are large.

In other words,

if both eigenvalues are small, there is little change in intensities along any direction => smooth region
if one eigenvalue is large and the other is small, there is change in intensities along one direction => edge
if both eigenvalues are high, then there is change in intensities along 2 perpendicular directions => corner

References:
https://stackoverflow.com/a/23172317/3337089
TODO:

Add proof showing positive definite matrices have positive eigenvalues
Add details about the curve trace - $\alpha$ * determinant


Difference between $\tt rand()$ and $\tt randn()$ in MATLAB rand() is for "Uniformly distributed pseudorandom numbers"

I generate two rand arrays, and then use the cross-correlation (xcorr2) and the normalized cross correlation (normxcorr2) as follows:
a = rand(512,512);
b = rand(512,512);
c = xcorr2(a,b);
figure;surf(c);shading flat;axis tight;
d = normxcorr2(a,b);
figure;surf(d);shading flat;axis tight;

The results are as follows:

in which there is strong correlation in the xcorr2().
randn() is for "Normally distributed pseudorandom numbers"

The same as the above:
a = randn(512,512);
b = randn(512,512);
c = xcorr2(a,b);
figure;surf(c);shading flat;axis tight;
d = normxcorr2(a,b);
figure;surf(d);shading flat;axis tight;

and the results are:

Then, my question is why there is apparent difference of xcorr2() for the two random data?

The reason for the difference is that rand() has an output range of $[0,1]$, and therefore you have a constant in both $a$ and $b$ in the first example. The other function randn() generates an output with zero mean. The constants are correlated. Try

a = -0.5 + rand(512,512);
b = -0.5 + rand(512,512);

or

a = rand(512,512);
b = rand(512,512);
a = detrend(a,'constant');
b = detrend(b,'constant');

and see if that is more like what you would expect...


How to Estimated the Noise Standard Deviation (STD - $ \sigma $) as a Function of Brightness from an Image? Hi I don't know where to start from.

I'm trying to implement a plot of noise with MATLAB.

I have found a plot of noise standard deviation sigma vs. brightness from a  images in the blog post - Profiling Sensor and Photon Noise.

But I think I need to some reference to implement profiling as above.

I also encountered the following reference - Noise Estimation from a Single Image.

How to make a plot of noise standard deviation sigma vs. brightness from a images?

I didn't read the article you referred to but I can try giving you some idea.

Run along the image for each pixel consider its 9 x 9 neighborhood.
For each pixel's neighbor hood calculate the STD and the Mean.

Use the Mean as the brightness value.
Plot the STD as a function of the mean value.

Some ideas for improvement:

Set the neighborhood size in an adaptive manner to only include data from the same brightness level.
For patches with the same value, aggregate the STD and average it.


Recoverable continuous function after average sampling I have a personal interest in this question, but I don't know if it is unfounded.
Suppose I have a continuous function x (t), which is sampled using an average sliding window, which performs a continuous uniform integration. The start of the scan is not fixed, that is, the signal could be scanned at any time. Is it possible to generate a function that is resistant to this procedure and can be recovered with an acceptable error from the samples obtained? What are the limitations in the generation of this function?
Thank you!

A rectangular prefilter is a sinc(x) in the frequency domain, thus quite far from the text-book brick-wall filter presented in standard discrete sampling.

Still, if your (constructed) continous waveform contains no frequency components at >= samplingrate/2 and you can compensate for the high frequency roll off, I would think that near perfect reconstruction is possible.

By multiplying the spectrum of your continous waveform with that of the fourier transformed rect (sinc), you should be able to see how far from classical generic sampling you are.

Incidentally, I think that both audio A/D (time domain) and camera sensors (spatial domain) really does something close to a sample-period integrator/sample-and-hold (quite rect-like) prior to sampling. Presumably because this maximize the signal-related energy that can easily be harvested above the noise floor?

-k


Sampling with an alternating impulse train The have the following question:

A signal $m(t)$ with bandwidth 500Hz is first multiplied by a signal
$g(t)$ where $\displaystyle g(t)=\sum_{k=-\infty}^{\infty}(-1)^k \delta(t-0.5*10^{-4}k)$. The signal is then passed through an ideal low
pass filter with bandwidth $1$kHz. The output of the low pass filter
would be?
The solution is given in the below image:

I am slightly confused with $G(f)$. As we know after sampling of a message signal $x(t)$ with frequency $f_m$ with the sampling rate being $f_s$, the output of the sampled signal $x_s(t)$ should have the following frequency components:
$$nf_s \pm f_m$$
Where $n=0, 1, 2, ....$
Then why according to the solution given there is nothing at $n=0$ i.e. $0 \times f_s \pm f_m=\pm f_m$, why there are no $\pm f_m$ components in the sampled signal or in other words why there in no component at $f=0$ for $G(f)$?
I was wondering if something is different due to the alternating impulse train.

We are not simply sampling the message in this case but also multiplying by a sampled tone that exists at $f_s/2$, (where $f_s$ is the sampling rate) and therefore we are frequency shifting the signal to this new location.

This is intuitively clear from the two plots above. In the first case, representative of a sampling process, we are multiplying our waveform by an impulse train in time. Such an impulse train in time is a impulse train in frequency, with each impulse separated by 1/T where T is the time spacing between the impulses in time. Multiplication in time is convolution in frequency, so the process replicates the spectrum of the waveform at each of the impulses shown (as you showed you understood from the sampling process). Given the impulse that exists at frequency location 0, there will be a copy of the waveform spectrum here that is unchanged from the original, assuming no aliasing occurs from the other replicas at the other impulse locations in frequency (convolving with an impulse at freq = 0 does nothing).

In the lower plot we are doing something similar, but in addition we are multiplying the waveform by a alternating +/-1 sequence. Note that such a sequence in time has no DC value! The plot is very clear if you consider the time domain waveform to be a sampled cosine waveform, with a frequency that is exactly half the sampling period T. So in this case the waveform spectrum will be shifted away from the DC location; to each of the impulses as shown.

As the OP has found, see this link that further details the mathematics: Expression to represent alternating pulse train?


What is the least JPG-compressible pattern? (camera shooting piece of cloth, scale/angle/lighting may vary) I am trying to design a cloth that, from the point of view of a camera, is very difficult to compress with JPG, resulting in big-size files (or leading to low image quality if file size is fixed).

It must work even if the cloth is far away from the camera, or rotated (let's say the scale can vary from 1x to 10x).

Noise is quite good (hard to compress), but it becomes grey when looking from far, becoming easy to compress. A good pattern would be kind of fractal, looking similar at all scales.
Foliage is better (leaves, tiny branches, small branches, big branches), but it uses too few colors.

Here is a first try:

I am sure there are more optimum patterns.
Maybe hexagon or triangle tessellations would perform better.

JPG uses the Y′ Cb Cr color space, I think Cb Cr can be generated in a similar way, but I guess it's better to not use uniformly the full scope of Y' (brightness) since camera will saturate the bright or dark areas (lighting is never perfect).

QUESTION: What is the optimum cloth pattern for this problem?

Noise is quite good (hard to compress), but it becomes grey when looking from far, becoming easy to compress. A good pattern would be kind of fractal, looking similar at all scales.

Well, there is fractal noise.  I think Brownian noise is fractal, looking the same as you zoom into it.  Wikipedia talks about adding Perlin noise to itself at different scales to produce fractal noise, which is maybe identical, I'm not sure:

I don't think this would be hard to compress, though.  Noise is hard for lossless compression, but JPEG is lossy, so it's just going to throw away the detail instead of struggling with it.  I'm not sure if it's possible to make something "hard for JPEG to compress" since it will just ignore anything that's too hard to compress at that quality level.

Something with hard edges at any scale would probably be better, like the infinite checkerboard plane:

Also something with lots of colors.  Maybe look at actual fractals instead of fractal noise.  Maybe a Mondrian fractal?  :)


literature related to Spectral resolution I  am interested in spectral resolution of spectral estimation methods.

Can you point me to some good literature, videos or books from which i can read and study?

I want to understand how to resolve  two closed spaced sinusoids:
$$
x_1(t) = \sin(\omega_1 t + \phi_1) \\
x_2(t) = \sin(\omega_2 t + \phi_2) \\
$$
where $\omega_1 \approx \omega_2$.

I would recommend
Array Signal Processing - Johnson, Dudgeon. This book covers classical spectral estimation, Minimum Variance Distortionless Response (MVDR), Linear Prediction, and subspace methods (e.g. MUSIC and ESPIRIT). It provides examples of the resolution capabilities between these methods. The are quite a number of good references in this book if you require additional detail.

MUSIC and ESPIRIT type algorithms are dicussed in most advanced DSP texts:
Advanced Digital Signal Processing, Proakis, Rader, Ling, Nikias
Discrete Random Signals and Statistical Signal Processing, Therrien

For parametric methods:
Spectral Analysis - A modern Perspective (Kay, Marple) Proc. IEEE Vol 6, No 11, 1981
Modern Spectral Estimation: Theory and Application, Steven Kay
Digital Spectral Analysis - S.L. Marple
SPECTRAL ANALYSIS OF SIGNALS - P Stoica and R Moses  PDF is available here

The books by Marple and Kay tend to focus on algorithms i.e. implementations of various MA, AR, and ARMA approaches (e.g. fast Lattice filter implementations) rather than the performance of a basic approach

Compressive Sensing / Sparse Reconstruction techniques can also be used. There are too many references here to list. There are a few books by these authors Michael Elad, Yonina Eldar, and Holger Rauhut. These techniques are often quoted as "having super resolution properties"

The algorithms are often hard to compare directly because they all have twiddle factors which affect their performance.  For example
 1. For MA, AR, ARMA, what order of model to use?
 2. For Compressive techniques - how closely to space the template vectors? What level of regularization to use? How many frequencies are there? How many samples do we have?


Signal spectrum when carrier width changes If I have spectrum $X_k$ of the signal evaluated at frequencies $k \delta f$ and I want to know what is the value of that same signal at some other frequencies , i.e., $k' \delta f'$.I know that $\delta f$ and $\delta f'$ are somehow reelated (e.g. $\delta f = 12 \delta f'$) how can I express spectrum of the signal with this new carrier spacing ($k' \delta f'$) having only info about $X_k$ or in other words what is the coefficient, $X_k' = coefficient \cdot X_k$.
Thanks in advance for any help.

If the signal is band-limited (at least above some noise level) below half the highest frequency in your initial frequency evaluation set, then you can resample the spectrum to any other frequency spacing by using periodic Sinc (Dirichlet) kernel interpolation.

If there is an integer relationship between the two carrier frequency bin spacings, then you can pre-calculate a polyphase interpolation table using the Sinc function (or a windowed approximation, if that meets your S/N requirements).  Otherwise you can convolve the entire spectrum vector against a computed periodic Sinc or Dirichlet kernel with the appropriate phase offset per computed "carrier" frequency.


HDR Imaging - Exposure Compensation Tried HDR Algorithms of OpenCV. Since the algorithm takes pretty much time. I tried few alternative steps as follows.

Took sample images from: HDR Images

1. Divided the image into four equal regions as Top Left and Right, Bottom Left and Right
2. Estimated the brighter region and the darker region from the Highly exposed image by taking average of the regions respectively
3. Replaced the brighter region with the average of Low and Mid exposed image.
4. Replaced the darker region with the average of High and Mid exposed image.
5. Then merged the images back into a single image.

 Can someone suggest how to compensate the difference in them
 Tried Exposure Compensator of OpenCV but does not get expected result.

 The result image looks like:

Here is one idea:
Merge the HDR images per-pixel using known thresholds for saturation etc. The final image will have pixel values above 256.

Apply a CLAHE (Contrast Limited Adaptive Histogram Equalization, its in OpenCV) to get the final image.

Often just doing the later will make your images look good.


Very simple question regarding $X(f)$ vs. $X(j\omega)$ Let's take a time domain function $x(t) = \cos( 2 \pi f_0 t) $. Its Fourier transform can be represented as

$$X(f) = \frac{1}{2} \left[ \delta(f - f_0) + \delta(f + f_0) \right]\tag{1}$$

as well as

$$X(j\omega) = \pi \left[ \delta(\omega - \omega_0) + \delta(\omega + \omega_0) \right]\tag{2}$$ which was initially wrong, corrected for clarity.

Considering $ \omega = 2 \pi f $, I can interchange between the two expressions. I can tell that the $j$ in the argument of the second expression can be accounted for by the $j$ in the RHS of the first one but I can't seem to get $(1)$ from $(2)$.

Starting with $(2)$:

\begin{align}
X(j\omega) &= \pi \left[ \delta(\omega - \omega_0) + \delta(\omega + \omega_0) \right]\\
X(j 2 \pi f) &= \pi \left[ \delta( 2 \pi f -  2 \pi f_0) + \delta( 2 \pi f +  2 \pi f_0) \right]
\end{align}

I can't see how to get $(1)$ from here. I know it is very simple but I can't get any good results by Googling it up. I'd appreciate if you do not downvote my post.

Your first equation is right but second equation is wrong:

I would proceed in the following way to see their relation. Given $x(t) = \cos(\omega_0 t)$ with $\omega_0 = 2 \pi f_0$, first see that:
$$x(t) = \cos(\omega_0 t) = \frac{1}{2}e^{j\omega_0 t} + \frac{1}{2}e^{-j\omega_0 t}$$
whose CTFT is:
$$X(\omega) = \pi [\delta(\omega - \omega_0) + \delta(\omega + \omega_0)] ..... (1)$$

which follows from the transform pair: $e^{j\omega_0 t} \leftrightarrow 2\pi \delta(\omega - \omega_0)$

Then use the following property of the impulse $\delta(a x) = \frac{1}{|a|}\delta(x)$ applied with $a=2\pi$ and $x=f$, and then (1) becomes:

\begin{align}
X(f) &= \pi [\delta(2\pi f - \omega_0) + \delta(2\pi f + \omega_0)]\\
X(f) &= \pi [\delta(2\pi (f - \omega_0/{2\pi})) + \delta(2\pi (f + \omega_0/{2\pi})]\\
X(f) &= \pi [\frac{1}{2\pi}\delta(f - f_0) + \frac{1}{2\pi}\delta(f + f_0)]\\
X(f) &= \frac{1}{2}\delta(f - f_0) + \frac{1}{2}\delta(f + f_0)\\
\end{align}
is the conclusion.

Note that the $j$ in the argument of CTFT $X(j\omega)$ is a notation that is equally valid in both forms as in $X(jω)$ or $X(ω)$, and it can be omitted for simplicity if there is no possibility of confusion, but otherwise remembered when it is necessary to underline the fact that the function is in general complex valued. Note that for the communications engineer the notation $X(f)$ is almost always is the choice. In either case arguments must be interpreted correctly based on the fundamental definition of the corresponding Fourier transform as given by:

$$X(jω) \triangleq \int_{-\infty}^{\infty}x(t)e^{−j\omega t}dt$$ , and
$$X(f) \triangleq \int_{-\infty}^{\infty}x(t)e^{−j2\pi f t}dt$$


Equalization knowing the Channel State Information I am trying to simulate a simple communication system using channel estimation. I was able to estimate the channel response using the Least Squares channel estimator, and it is working properly (evaluated by the mean square error between the channel and its estimation). But once I have this knowledge, I need to equalize the signal at the receiver. I have tried to apply the zero-forcing method, but it's not providing a good performance.

My questions are:

What is the proper way of doing this equalization?
In a system where two transmitters send information at the same time and in the same frequency, I can estimate both channels at the receiver, but how can I equalize the signal once I know both channels?
Is there is any reference that anyone could indicate to read more about this topic and help answer these questions? It wold be enough for now.

Thanks in advance.

To 1.
From the channel impulse response you can build the channel matrix and then compute the ZF equalizer coefficients using the pseudo inverse of the channel matrix.

Once you have the ZF coefficients you can convolve them with the received signal.

To 2.
I guess you are referring here to a PNC system. The best solution is to use OFDM to solve this problem. In this case, you can use Maximum Likelihood detection per subcarrier by checking which combination minimises the Euclidean distance.

To 3.
You can find the equations for ZF equalizer in the book by J. Kurzweil, An Introduction to Digital Communications, 1st Edition.

BTW, the performance of the time-domain ZF equalizer is not very good in frequency selective channels. You are better off with the decision feedback
equalizer.


Strange Response of a fir1 Filter in Matlab i have been tinkering with some fir filters to bandpass-filter a signal. Even though the filter looks pretty descent in frequency domain, i receive a quite strange "filtered" signal in time-domain. I think, i'm oversseing something simple... Maybe one of you knows what i am doing wrong.

First of all, here's the code:

Fs          = 200; % hz
nyquistFreq = Fs/2;
frange      = [45 55];
order       = round( 30*Fs/frange(1) );
filtKern    = fir1(order,frange/nyquistFreq);
filtSignal  = filter(filtKern,1,signal);

I plotted the data stored in signal and filtSignal into one plot and the filter-shape is received at the beginning of the filtered signal:

Here's a picture of the filter in time-domain, in frequency-domain and in frequncy-domain with logarithmic scale. According to the plot in frequency domain, the cut-off frequencies fit pretty much what i requested...

Thank you so much, i appreciate every input!

Since you are seeing your filter coefficients and much larger than your signal, this is indicative that you have a very large "impulse" at the start of your signal. The coefficients of your filter is the impulse response for the filter, so that is exactly what you are seeing: the response to an impulse.

Review the start of your time domain data for a very large offset in the first sample or near the start of the file.


How to design a digital filter with a frequency dependent phase delay? I want to design a filter with a custom phase delay related to frequency. As frequency increases, phase delay should increase.

The time delay as a function of frequency can be expressed as:
$t_d = \frac{L}{C_{ph}} - \frac{1}{f}$

$L$ and $C_{ph}$ are constants of $0.0017$ and $2628$ respectively.

The range of $f$ is $500\mathrm{kHz}$ to $1 \mathrm{MHz}$, the sampling frequency, $f_s$ is $78.39\mathrm{MHz}$.

Thus, when the frequency is $500 \mathrm{kHz}$, the delay should be $1.37\mathrm{\mu s}$ or $107$ samples.

For the pass band, $f$, I have used a filter response of $e^{-i2\pi fd}$, where $d$ is the delay in samples required.

I have designed stop bands with a filter response of zero between $0 \mathrm{Hz}-100\mathrm{kHz}$ and $1.4 \mathrm{MHz}-f_s$.

In MATLAB my code looks like this:

n = 50;
fs = double(fs);
L = double(L);

Cph = 2628;

f1 = linspace(500e3/fs,1e6/fs,100);
f1d = L/Cph - (1./(f1*fs));
%Our array is backwards though innit
f1d = f1d * -1;
f1dz = f1d * fs;
h1 = exp(-1i*2*pi*f1.*f1dz);

fstop1 = linspace(0,100e3/fs, 10);
hstop1 = zeros(size(fstop1));

fstop2 = linspace(1.4e6/fs, 1, 10);
hstop2 = zeros(size(fstop2));

d=fdesign.arbmagnphase('n,b,f,h', n, 3, fstop1, hstop1, f1, h1, fstop2, hstop2);
%d=fdesign.arbmagnphase('n,b,f,h', n, 1, f1, h1);
D = design(d,'equiripple');

fvtool(D,'Analysis','phasedelay');

This is what I get:

Markers are at 500k and 1M.

What am I doing wrong?

I notice two problems in your code:

at $f=500\,\text{kHz}$ you require a phase delay of $107$ samples, but your filter length is only $50$ taps; you can't have a delay that's greater than the filter length of an FIR filter.
it looks like you normalize all frequencies by the sampling frequency; however, in Matlab a normalized frequency $f=1$ corresponds to the Nyquist frequency, i.e., $f_s/2$, not to $f_s$.


Linearity of the given system I am given the following system and I am checking the additive property:

$$y(t)=x(e^t)$$

where $y(t)$ is the output and $x(t)$ is the input given to the system.

Now this is what I did so far:

\begin{align}
x_1(t)&\rightarrow \text{system} \rightarrow& y_1(t)=&x_1(e^t)\\
x_2(t)&\rightarrow \text{system} \rightarrow& y_2(t)=&x_2(e^t)\\
\end{align}

$$y_1(t) + y_2(t) = x_1(e^t) + x_2(e^t)$$
How will system respond to this input?

$$x_1(t) + x_2(t)\rightarrow system \rightarrow = $$

What I think is that this is what will happen:

$$x_1(t) + x_2(t)\rightarrow system \rightarrow = x_1(e^t)+x_2(e^t)$$

Am I right?

Additivity requires a little more than a direct addition of some $x_1$ to some $x_2$. It should involve any linear combination of  $x_1$ and $x_2$, $a_1 x_1 + a_2 x_2$ where $a_1$ and $a_2$ are "kind of scalars" (this could make us dive into complicated stuff, like in  Does scaling property imply superposition?)

In your case, however, linearity works, because you are touching the "time base" ($t$ vs. $e^t$) without touching the amplitudes. Your operation is a sort of time warping.

Accordingly:

$$\mathcal{S}(a_1 x_1(t) + a_2 x_2(t)) = a_1 x_1(e^t) + a_2 x_2(e^t) = \mathcal{S}(a_1 x_1(t)) + \mathcal{S}(a_2 x_2(t))  $$


Minimum Sampling frequency for this signal to avoid aliasing I would like to compute the minimum sampling frequency for the signal $x(t)$ given below:

$x(t)= sinc ^3(300t) * sinc ^5(600t)$ and  $ *$ stands for convolution.

where the sinc function is defined as $sinc (w t) = sin (\pi wt)/(\pi wt)$.

I have tried like this  : $sinc^3(300t)$ will correspond convolution of 3 rectangular functions in frequency domain which results in a frequency limited function to 3 times 150 Hz ie 450 Hz. And the other function corresponding to maximum of frequency 5 times 300 Hz which is 1500 Hz. Now finally convolution operator get transformed to multiplication in frequency domain , is the maximum frequency content limited to 450 Hz and minimum sampling frequency 900 Hz?

It would be helpful if someone shows me the derivation of fourier transform of sinc function for this purpose.

There's several different ways you can go about solving this problem just in the time domain. Use trig identities or Eulers formula to reduce the equation. Both are simple enough and just require some time to work through. Feel free to use WolframAlpha to check your work. I find it sometimes helpful for tedious algebra.

Yes, the Fourier transform of a sinc function is a rectangular function, so you will have what you've described as convolved rectangular functions. This is more work than is necessary as after you've simplified the time-domain expression you will be able to find the maximum frequency, and subsequently the required sampling rate, for that signal.


Bandwidth and Bit rate I'm kind of confused about digital transmission.

Is the bandwidth occupied by some digital signal the Fourier spectrum of the bit (pulse) format,or the one determined by the bit rate?

Bandwidth is of course occupied by the signal itself. The modulation design of the signal determines the bandwidth. That's why pulse shaping techniques are adopted  by many different kind of modulation techniques.

However, our bandwidth we occupy has a DIRECT relation with our maximum bit rate. In Shannon's theorem:

$$
C = B * log_2(1 + S/N)
$$

where

C = channel capacity (bits/s),
B = channel bandwidth (Hz),
S = signal strength (watts),
N = noise power (watts).

This theorem tells us that we CANNOT transmit more bits per second through that channel in good condition.

Hence, our signal determines our bandwidth, and our bandwidth determines our maximum achievable bit rate.

It depends on some number of things like Inter Symbol Interference and Nyquist Theorem. It's relatively a long topic to explain each element. Please take a look at the terms I have mentioned (Nyquist and Inter-Symbol Interference) and if there is a spot if you have a hard time to understand, please don't hesitate to ask.

As always, Have a nice day.


Are time series data always contain noise? I wonder if every time series data should contain noise or not. For example I am taking the price of a ticker, say Yahoo, every hour and noting the values. Does this data contains noise or not?

Periodically taken ticker data is very likely not a set of samples of a band-limited function, thus can contain aliasing artifacts which can be considered noise.  This is possibly a much larger contribution than quantization and time jitter noise.

Absence of noise requires the existence of well defined signal being measured.  Proving that such a signal really exists is a separate and likely far more difficult problem for some types of time series.


On the application of chaos in signal processing Paper: A novel channel equalizer for chaotic digital communication system (DOI, RS link).

I am having problem in understanding their claim which is that a nonlinear system, considered here to be a chaotic system given by Eq (1)
$x[n] = f(x[n-1])$ is represented as a symbolic system $x[n] = 0.b_nb_{n+1}...$ where the $b$'s are {0,1} and so a time series in terms of the binary expansion of a number $x[n]$ describes the dynamics of the system. This is explained in Section 2. Then they say that a linear filter of the form Eq(5)

$x[n] = \sum b[n-k]h[k]$

where the $b$'s are the information bit sequence and the output $x[n]$, is chaotic. So, this Equation according to them is a way of linearizing the nonlinear dynamics of Eq(1). If this is a linear representation, then how come it is chaotic?

I am unable to understand how the bits-- $b$'s that they are talking about is generated and how the binary bit sequence can be used as the input
to the filter of Eq. (5) to get a chaotic data sequence i.e. why
is this representation chaotic? One of the ways to produce symbolic dynamics of a trajectory of a map is by thresholding explained as:

b[n] = 1, if x[n] >= threshold, otherwise b[n] = 0

Is this how the $b$'s are generated from the same map or these bits are coming from some other source?

A typical communications system is composed of the following parts:

Message (Information Source) --> Encoding --> Modulation --> Channel --> Demodulation --> Decoding --> Recovered Message

A "chaotic modulator" is just another way of doing the Modulation building block. This is depicted in Figure 1 of the paper that is cited in the question.

The input to the Modulator is a binary sequence and the output of it is a digital (or directly analog, depending on the system) signal.

Therefore, to clarify (rather than answer) the first issue: As per Figure 1, the bs are essentially your Message (or Message after it has been encoded for whatever purpose (e.g. error correction).

A typical Modulator's job is to map binary symbols to carrier symbols with the carrier signal being one that is suitable for transmission effectively through the Channel. In a typical system, we can literally see this mapping in the Constellation Diagram.

In the case of a chaotic modulator though, binary symbols are mapped to different states of the chaotic modulator. This is mentioned in the introduction of the paper as:

  Rather than using structured signals, such as rectangular  pulses or
  sinusoids, to denote ‘0’s and ‘l’s, these communications systems embed
  the information in the time evolution, or dynamics, of the transmitted
  signal.

A simpler way to think about this is to take some difference equation (such as f) and produce a few values with some initial condition x0 and a few values for some initial condition x1. If the two time series are distinct enough, then perhaps it is possible to recover which of the two x0,x1 might have generated this time series. In this case, of going from time series to some initial condition, we are looking at some sort of Demodulator.

Now, chaotic of not, a difference equation provides a model. A model imposes a set of constraints and some sort of structure to the signal. This is what is exploited here (please keep in mind that we are talking about a pseudo-chaotic system). So, what if I take look at the difference equation and realise that I can express it as a typical digital filter?

Again, to clarify, rather than answer, this is what Drake & Williams did with "Pseudo-chaos for direct sequence spread spectrum communications". From the paper linked in the question:

  More  importantly, Drake  and Williams  [7] have  shown  that the
  sawtooth  map's  output is equivalent  to the response of a linear
  filter to a binary sequence.

Having clarified this, it is also important to highlight Figure 2 (and 3 to some extent) from the original paper which shows that by inversing the filtering operation it is effectively possible to do Modulation and Demodulation.

This, hopefully, clarifies the "representation" missunderstanding.

Finally, I would like to stress the pseudo-chaotic nature of the generated signal which comes as a result of the modifications that are also explained in section 2, from "In practice this is not a possible because this filter..." (sic) onwards.

I would also like to recommend that you take a brief look at Pseudo-random number generators, the Linear congruential generator and finally the linear feedback shift register (LSFR) for examples of seemingly linear "representations" that still produce rather complex sequences and then taking a brief look at the Dyadic Mapping where it is also mentioned that:

  The name bit shift map arises because, if the value of an iterate is
  written in binary notation, the next iterate is obtained by shifting
  the binary point one bit to the right, and if the bit to the left of
  the new binary point is a "one", replacing it with a zero.

Hope this helps.

EDIT:

In view of additional comment:

  If the chaotic map is denoted by the state variable $x$, then these
  $b$ are the bits. Further, are the $b$′s obtained from the binary
  expansion of some numeric value of the iterate of chaotic map $x$?

Here is an alternative way to go about this:

Let's say that this is not a modulator. Ultimately, $x$ is the output from the mapping. We start with some initial value $x_0$ and produce $x_1$ which we send back into the map and get back $x_2$ and again through the map to get $x_3$ and back through it once more to get a $x_4$ and so on to get $x_n$. Since we are operating on a computer with numbers of finite wordlength, $x$ contains DIGITAL values (in the decimal system) with a well defined minimum and maximum. So, if $x$ is a byte, then the min value is 0 and the maximum value is 255.

Each one of these $x_n$ decimal values has a binary representation as well. Therefore:

$x_0$ is equivalent to some binary number $0100101$

$x_1$ is equivalent to some binary number $1001010$

$x_2$ is equivalent to some binary number $0010100$

And eventually...

$x_n$ will be equivalent to some $b_m b_{m+1} b_{m+2} ...$ where some $b \in {0,1}$

Therefore, it is possible to have a binary string of digits (0010001101...) for a pseudo-chaotic sequence $x$.

Let us now ask the reverse question: What if all we have is the binary sequence. Can we go from the binary sequence to the output it would produce?

To do this, you would have to establish a relationship between the $b$ (binary sequence) and $x$ (decimal sequence / digital signal). It turns out that this is possible via the work of reference 7 as mentioned above, through $h$.

To recap and put it in the context of communications:

$b$ is your Message. If the message was "ABC", its ASC representation would be 65,66,67 and the binary representation of that would be "010000010100001001000011". This is your $b$.

$x$ is the output from the Modulator with pseudo-chaotic dynamics.

$h$ is now the tool by which the map is applied to the data.

(Please see figures 2,3 in the original paper).

  How can the convolution of $b$ with impulse give rise to a chaotic
  map, $x$?

Again...The input to the modulator is a binary sequence ($b$). Each $b$ is either $0$ or $1$. The output from the modulator is a digital signal ($x$) (decimal values of some wordlength).

If we take a look at:

$x[n] = \sum_{k=-\infty}^{\infty}{b[n-k] \times h[k]}$

We realise that:

$x$ is a weighted sum of $b$.
To produce one $x_n$, we need to sum many bits from $b$. How many? $k$ many.
Obviously, $k$ cannot remain infinity in a practical system. That's where truncation of $h$ comes in.

I hope this helps.


Seeing spurious peaks in acquired vibration data; peak frequency changes with sampling rate I am acquiring vibration data from an accelerometer using a DAQ. Due to some reason I am seeing spurious peaks/impacts in the FFT spectra at specific frequency and its harmonics all the way till the Nyquist frequency. This frequency changes with the sampling rate as follows:

Sampling rate = 48 KHz, peaks at 100 Hz and its harmonics
Sampling rate = 8 KHz, peaks at 83.33 Hz and its harmonics
Sampling rate = 24 KHz, peaks at 93.75 Hz and its harmonics

These peaks don't seem to correspond to any real mechanical/vibration source. Initially I suspected aliasing, but it seems the DAQ has anti-aliasing built-in. I also tried doing anti-aliasing in software using python/scipy.signal.decimate(), but it didn't help.
I am attaching below couple of FFT spectra plots for the 24 KHz case.
Can you please advise what could be possible reason for such spurious FFT peaks?

Can you please advise what could be possible reason for such spurious FFT peaks?

Most likely it's your measurement setup. Good  measurements are hard to make and require careful setup, characterization, and execution.

Determine the noise floor (spectral shape and level) of your data acquisition system, just short the input and grab few spectra. If that looks good
Determine the noise floor (spectral shape and level) of you sensor: attach the sensor to something heavy that's not moving and grab few spectra. If that looks good
Build a "calibration" setup. Attach your sensor to a well characterized, calibrated or programmable vibration source (a shaker, or maybe a electro mechanical driver) and measure different frequencies at different levels and compare extend the expected results. Determine the frequency and dynamic limits of your setup.

In each of the steps above set targets, requirements, and expectations BEFORE you make the measurement. These should come from the data sheets, general physical considerations and the specifics of your application.
Don't even look at any real measurement data until you have proven that your setup does indeed measure vibration with sufficient accuracy and signal to noise ratio over the required frequency range and dynamic range. If that all checks out and you still see unexpected data you can ask again.