Why can't a causal digital filter have an infinitely sharp transition between the passband and the stopband? In DSP book by Proakis and as well as in this pdf, it is mentioned that practical causal digital filters cannot have an infinitely sharp transition from Pass-band to Stop-band. Why is it so? Can you please provide a detailed explanation (with proof)?
Edit 1: In the book in was just mentioned that is a consequence of the Gibbs phenomenon, which results from the truncation of h(n) to achieve causality. I didn't understand how.

I don't have a concrete proof for this one. However, I can tell you this...
Consider a perfect low pass filter. The time domain representation is a sinc. And for any system to have a sharp transition band, a base signal has to be multiplied with a rectangular waveform in the frequency domain. Which implies that, the time domain signal of the same has to be convoluted with a sinc in time domain. We know that sinc function is a non causal signal.
Hence proved.


Is sift a good way to extract features from an image? I am trying to extract features from an image, but I have failed to get the points that I want to extract, and hence my image fails to match with the template.

Are there any checkpoints that I have to go through before applying SIFT to my image so that I get good results? For example the template image is;

targeted image;

What you are supposed to do when matching a template to an image using sift is to run sift against your template and then look for those sift features in that arrangement in your scene.

Rule of thumb: Compare like to like.

Sift(Template) Contained Within Sift(Image)

You cannot tune Sift to extract the "features you want" Sift *uses invariant measures to find features that it believes are most salient. (*Chris A)

If you want it to find particular features from your template make smaller templates of those features and search for them.

Just a suggestion. I do not know exactly what type of template you are trying to fit to what type of scene or image.


Linear phase Notch filter / Band reject filter implementation in C++ I'd like to implement a digital notch filter, which has linear phase, in C++ or C#.
Example : a notch filter that removes 400hz frequency on a .wav file (16 bit, 44.1 Khz, stereo).
I'm a bit lost with many examples, I don't achieve to implement them.
Maybe some of you has a working example?

Thank you in advance.

You can implement a simple recursive filter to process blocks of samples like this:

void filter(const int *x, int *y, int n)
{
    static float x_2 = 0.0f;                    // delayed x, y samples
    static float x_1 = 0.0f;
    static float y_2 = 0.0f;
    static float y_1 = 0.0f;

    for (i = 0; i < n; ++i)
    {
        y[i] = a0 * x[i] + a1 * x_1 + a2 * x_2  // IIR difference equation
                         + b1 * y_1 + b2 * y_2;
        x_2 = x_1;                              // shift delayed x, y samples
        x_1 = x[i];
        y_2 = y_1;
        y_1 = y[i];
    }
}

Notch filter coefficients a0, a1, a2, b1, b2 can be calculated using eq 19-8 on this page: http://dspguide.com/ch19/3.htm


Getting frequency content at different times from discrete wavelet transform coeffs After being away from DSP for a long time, I am trying to familiarize myself with wavelet transform. Here is what I (think) have understood so far:

Wavelet transform provides you high time resolution at higher frequencies and high frequency resolution at lower frequencies.
DWT can be calculated by using QMF pair and subsampling. When used recursively, filter pair increases frequency resolution and subsampling decreases time resolution.
Result of level 1 DWT [cA, cD] = dwt(signal, Lo_D, Hi_D) essentially gives low frequency and high frequency splits of the signal called approximation and details.

If my understanding is correct, in this algorithm, the end result is essentially in time domain. If I want to know what frequencies are present, I have to take an fft of the coefficients. Is that correct? If so, is there a way I can get the frequency domain representation without extra step of taking fft?

Question 2: If scalogram is the answer to question 1, how is it generated using results of repeated use of dwt?
I want to avoid using ready MATLAB functions for better understanding except maybe dwt

Thanks

Question 1: I think you want to investigate centfreq &  scal2freq (Scale to frequency):

  FREQ = centfrq('wname') returns the center frequency in hertz of the
  wavelet function

  F = scal2frq(A,'wname',DELTA) returns the pseudo-frequencies
  corresponding to the scales given by A and the wavelet function
  'wname' (see wavefun for more information) and the sampling period
  DELTA.

The intuition that I use is to think about your question is: The wavelet transform is a hierarchical (tree-like) decomposition of an input signal. I would not say it is in the "time-domain" as is commonly understood. I think of the DWT output via the duality: Wavelets <---> Filterbanks.

Imagine a set of 8 bandpass filters tuned to filter different ranges from low to high on the frequency spectrum. You know which frequency ranges correspond to which filters, so you know how much signal lives in each of the 8 frequency ranges based on the bandpass filter output. That's a basic frequency analysis. You wouldn't need to compute an FFT at this point. And with the successive QMF filtering, then decimation of your signal at each level, your DWT coefficients may only be one-sample (or have a very coarse resolution), so an FFT wouldn't make any sense.

There's a good explanation @ https://math.stackexchange.com/questions/28581/which-time-frequency-coefficients-does-the-wavelet-transform-compute

Question 2: A scalogram (Scalogram (and related nomenclatures) for DWT? ) is good for visualizing your transformed signal, but I think it would be overkill just to compute which frequency range your wavelet coefficients correspond to.


DC removal from IQ samples I am wondering what is the proper way to remove DC component from PSK modulated signals. Should I calculate mean value for I and Q channel separately and correct then samples in each of the channel or calculate one mean value for I and Q channel as one channel with interleaved samples (IQIQIQIQ) and then correct samples. What do you think?

"What do you think" is the wrong question here: It's not a matter of opinion. Your second approach is mathematically simply not calculating the mean of a series of complex numbers. The mean of a complex sequence is w.l.o.g. complex, the mean of your IQIQIQIQ is real, so it can't possibly be the mean/DC component.


How can I determine the order of filter? This is my first study about signal analysing. I'm very confused about filter order. My problem is how can I know whether its 12-th order, or 2nd order like the book says so? I already knew that slope has 2 choises whether its in n multiples of -6n dB/oct or -20n dB/dec. If the slope is -20 dB/dec then the second order should have -40 dB/dec. For example, An n= 4th order is 24dB/octave slope. But, I saw this example below and pretty confuse about determine the slope of those two line and how can the first line is 12-th order. Its not like the slope reaches -240 dB/dec (-20 x 12) and the second line seems not like -40 dB/dec (-20 x 2). Is there anyone could explain this to me? Thank you.

Note that the given slopes of $6n\textrm{ dB}/\textrm{octave}$ and $20n\textrm{ dB}/\textrm{dec}$, respectively, are asymptotic slopes for large frequencies. So you have to extrapolate the lines that are approximated by the magnitude responses at large frequencies.

For the given figure you could do the following. For the second order filter, look at the attenuations at $f=10^2\textrm{Hz}$ and at the next vertical line, which corresponds to $f=2\cdot 10^2\textrm{Hz}$. The ratio between those two frequencies is an octave, and the difference in attenuation is approximately $-22-(-35)\textrm{dB}=13\textrm{dB}$ (should be $12\textrm{dB}$).

For the $12^{th}$order system you could look at the attenuation values at $30\textrm{Hz}$ and $40\textrm{Hz}$, taking into account that the curve flattens towards its knee, and you need to extrapolate the straight line. With this in mind you would arrive at approximately $-31\textrm{dB}$ at $f=40\textrm{Hz}$, and at $-2\textrm{dB}$ at $f=30\textrm{Hz}$, which corresponds to a difference of $29\textrm{dB}$. Solving

$$\left(\frac{4}{3}\right)^x=2\tag{1}$$

for $x$ gives $x\approx 2.41$, i.e., you have to multiply $29\textrm{dB}$ with $2.41$ to obtain the attenuation per octave: $29\cdot 2.41=69.9\textrm{dB}$ (should be $72\textrm{dB}$).


m-ary antipodal signaling I know that the error probability for binary antipodal signaling is $Q(\sqrt{2E_b/N_0})$ with $Q$ the tail function.

I was wondering if when used in a m-ary scheme how the probability error would change. What I mean by that would be a signaling where we still have only 2 antipodal codewords, but in a ${\rm I\!R^m}$ dimension domain, i.e. $c_1$ = +$\sqrt{E_b/m}$(1,..1)
$c_2$ = -$\sqrt{E_b/m}$(1,..1)

Thank you for reading, and have a great day.

Let the received vector $\vec{y}$ be given as $$\vec{y} = \vec{x}+ \vec{w}$$ where $\vec{x}$ can be any of $c_1$ or $c_2$. The ML detection rule for deciding $c_1$ is given as $$ \frac{pdf_Y(y;c_1)}{pdf_Y(y;c_2)} > 1\tag{1}$$  $$pdf_Y(y;c_1) = (2\pi)^{-M/2}(\sigma^{-M})e^{-\frac{\Vert \vec{y} - \vec{c_1} \Vert^2}{2\sigma^2}}\tag{2}$$ similary $pdf_Y(y;c_2)$ is $$pdf_Y(y;c_2) = (2\pi)^{-M/2}(\sigma^{-M})e^{-\frac{\Vert \vec{y} - \vec{c_2} \Vert^2}{2\sigma^2}}\tag{3}$$ using $c_1 = -c_2$ and (2), (3) in (1) taking log of both sides and simplifying we get the detection rule for deciding $c_1$ as $$ \vec{y}^Tc_1 > 0$$.

Now an error will happen when  $c_2$ gets transmitted but $ \vec{y}^Tc_1 > 0$.

Therefore we have $$ (c_2 + \vec{w})^Tc_1 > 0$$ Now $z= \vec{w}^Tc_1$ is a gaussian random variable with mean 0 and variance $\frac{E_b}{M}M\sigma^2 = E_b\sigma^2$. Therefore we have $z > \Vert c_1\Vert^2$. Now normalizing $z$ with standard deviation of $z$ (to convert it standard normal with standard deviation 1 and mean 0) we have $$\frac{z}{\sqrt(E_b\sigma^2)} > \frac{\Vert c_1\Vert^2}{\sqrt(E_b\sigma^2)}$$ Now $\Vert c_1\Vert^2 = E_b$ therefore we have probability of error as $$Q(\sqrt\frac{E_b}{\sigma^2})$$ since $\sigma^2 = \frac{N_o}{2}$ (per unit bandwidth). Therefore the error probability is $$Q(\sqrt\frac{2E_b}{N_o})$$

The assumptions are that noise in $\vec{w}$ are independent and identically distributed.

For the simple case of binary antipodal signal we have an energy per bit of $E_b$ and noise in 1 dimension. Now for the extended case we have energy per bit of $E_b/M$. Hence the expression above

Note: If the noise along each dimension is not IID, then we would not have this same expression even though eucledian distance between the signal points remains the same. But you can use the approach I presented here to compute the bit error rate.


Noise removal from audio using FFT/Spectral editing I have a noisy recording with some letters spoken out loud. I've noticed that the noise is easily recognizable when looking at a spectrogram of the audio. Is there anyway to set a power/frequency threshold for the spectrogram, and muting the parts which do not reach that threshold? (The noise doesn't have enough power in higher frequencies, as you can see in the spectrogram)

There exist solutions that are coarsely based on the idea that you propose: Threshold the spectrogram amplitudes to reduce/remove the noise part of the signal. Usually, this will bring about artifacts as e.g. the so-called musical noise and the application determines whether this is acceptable or not. In general, the more the noise is reduced the more negative effects the process has on the quality of the desired signal (i.e. the speech in this case).

As a first attempt I would suggest you to try the Noise Removal algorithm in Audacity (http://web.audacityteam.org/). Probably this will give satisfactory results and is sure worth a try as the software is free. If you're okay with spending money to most probably achieve better results, have a look at RX by iZotope (https://www.izotope.com/en/products/audio-repair/rx/). However, if you know the basics of digital signal processing and want to experiment yourself, a first thresholding algorithm will not be very complex in programming environments like e.g. MATLAB.

If you're interested in the background / signal processing methods that these algorithms are based on I would recommend the books by Vaseghi (http://www.amazon.com/Advanced-Digital-Signal-Processing-Reduction/dp/0470754060/ref=sr_1_1?ie=UTF8&qid=1432211838&sr=8-1&keywords=vaseghi) or Godsill and Rayner (https://www-sigproc.eng.cam.ac.uk/Main/SJGSpringer), the second of which can be downloaded for free from the internet.


How do I average frequency spectra? How did letter d.) become letter e.)? How do we perform this averaging process? Its clear the averaging the signal causes some noise to be reduced.

If you are in a position to acquire $N$ snapshots of data, $s_n(t)$, for $n=0\ldots N-1$ then you can just do:

$$
S_n(k) = |\mathrm{fft}(s_n(t))|^2\\
\hat{S}(k) =  \frac{1}{N} \sum_{n=0}^{N-1} S_n(k)
$$

and a certain amount of averaging will reduce the apparent noise on the spectrum.


Destructive interference on the Autocorrelation of the time signal of periodic CPFSK-signals. An FSK-signal with a frequency-shift $\Delta F$, a symbol-stream $x$, $x(n) \in {\{-1,1}\} $ and symbol duration $T$ has the complex envelope:

$$f[x](t)=\exp\left(2 \pi i \Delta F T \left(\sum_{l=0}^{j-1} x(l)\right)+\left(\frac{t}{T}-j\right)x(j)\right)\quad\text{for}\quad jT\leq t\leq (j+1)T$$

I am examining periodical symbol-streams for FSK-signals and I found a peculiar dependency between the strength of the periodicity and the modulation index,which I am very unsure of:

The phase of the samples at the start of a symbol is

$$2\pi i \Delta F T \left(\sum_{l=0}^{j-1} x(l)\right)=:\varphi(j,T).$$

So the phase deviation between the samples at the beginning of two consecutive symbols is given by

$$|(\varphi(j,T) -\varphi((j-1),T)|=\pi \eta, \tag{1}$$

where $\eta=2 \Delta F T$ is the modulation index of the FSK signal.

Due to $(1)$ for integer modulation indices $\eta$ the phase at the beginning of a symbol stays constant. Accordingly, a complex autocorrelation of a periodic signal shows a peak at right at the length of a period of the signal.

But, if the modulation index is not integer, the starting phase of the symbols shifts throughout the signal, i.e. $2\pi \nmid \pi \eta $.
A complex autocorrelation of same period signals, shifted in phase, does not show the same peaks and therefore does not exhibit periodicity in the same magnitude.

Is my suspicion  really sound for continuous phase frequency shift keying signals?
Does this hold for MSK-signals too ?

After thinking for a while i came up with my own proof. I hope it is helpful to you and correct!

Definition: Frame-Structure

The frame structure used in the following is one of one-periodic signals, with a constant header and a uniform distributed payload. The ratio between payload and total frame is given as $$ p_{\%} = \frac{|payload|}{|frame|} $$
the ratio of samples in the header vs samples per frame.

Example:
$$ |HHH|012|HHH|223|HHH|978|...$$
,here: $p_{\%}=0.5$ .

Proof:

The complex envelope of an FSK-signal is

\begin{array}{cc}
x\left(t\right)=e^{i\left[2\pi\Delta FT\left(\varphi\left(j,T\right)+\left(\frac{t}{T}-j\right)d\left(j\right)\right)\right]} & ,\text{}jT\leq t\leq\left(j+1\right)T\end{array}

with

$$\varphi\left(j,T\right):=2\pi\Delta FT\sum_{l=0}^{j-1}d\left(l\right)+\varphi_{0}$$

and the modulation index $$\eta:=2\Delta FT.$$

We introduce the model by having a look at the ACF of truely periodic, i.e. $p_{\%}\left(x\right)=1 $, discrete signals. In a second step we generalize to partly periodic signals with $p_{\%}\left(x\right)<1$. Let $x$ be a digital FSK signal of length N+1, that is, $x$ is a sampled version of the complex envelope above. Then the ACF of x is defined by

\begin{array}{ccc}
\rho\left[x\right]\left(\tau\right) & = & \sum_{k=0}^{N-\tau}x\left(k\right)x^{*}\left(k+\tau\right)\\
 & = & \sum_{k=0}^{N-\tau}e^{i\left[\pi\eta\left(\left(\frac{k}{T}-j_{k}\right)d\left(j_{k}\right)+\varphi\left(j_{k},T\right)\right)\right]}e^{-i\left[\pi\eta\left(\left(\frac{k+\tau}{T}-\left(j_{k+\tau}\right)\right)d\left(j_{k+\tau}\right)+\varphi\left(j_{k+\tau},T\right)\right)\right]}
\end{array}

Let $x$ be $\tau_{0}$-periodic, consequently $\tau_{0}$ is an integer-multiple of the symbol duration T measured in samples. Given our data model, we can expect to see a distinct peak in the ACF feature of such a signal. Let $j_{k}:=\left\lfloor \frac{k}{T}\right\rfloor$  be the index of the symbol being modulated at point of time $k\in\mathbb{N}$. While the signal $x$ is $\tau_{0}$-periodic, the series of symbols is $\frac{\tau_{0}}{T}$-periodic, i.e. $j_{k+\tau_{0}}=\left\lfloor \frac{k+\tau_{0}}{T}\right\rfloor =j_{k}+\frac{\tau_{0}}{T}$. Therefore the symbols of a signal are repeated periodically, i.e. $d\left(j_{k+\tau_{0}}\right)=d\left(j_{k}+\frac{\tau_{0}}{T}\right)=d\left(j_{k}\right)$. Hence we can simplify the ACF of x at lag $\tau_{0}$ to

\begin{array}{cc}
\rho\left[x\right]\left(\tau_{0}\right)\\
= & \sum_{k=0}^{N-\tau_{0}}e^{i\left[\pi\eta\left(\left(\frac{k}{T}-j_{k}-\frac{k+\tau_{0}}{T}+j_{k}+\frac{\tau_{0}}{T}\right)d\left(j_{k}\right)\right)\right]}e^{i\left[\pi\eta\left(\sum_{l=0}^{j_{k}-1}d\left(l\right)-\sum_{l=0}^{j_{k}+\frac{\tau_{0}}{T}-1}d\left(l\right)\right)\right]}\\
= & \sum_{k=0}^{N-\tau_{0}}e^{i\left[\pi\eta\left(\left(-\frac{\tau_{0}}{T}+\frac{\tau_{0}}{T}\right)d\left(j_{k}\right)\right)\right]}e^{\pi\eta i\left(-\sum_{l=j_{k}}^{j_{k}+\frac{\tau_{0}}{T}-1}d\left(l\right)\right)}\\
= & \sum_{k=0}^{N-\tau_{0}}e^{\pi\eta i\left(-\sum_{l=j_{k}}^{j_{k}+\frac{\tau_{0}}{T}-1}d\left(l\right)\right)}.
\end{array}

The value $\rho\left[x\right]\left(\tau_{0}\right)$ is decisively influenced by the term $p_{m}:=-\sum_{l=j_{k}}^{j_{k}+\frac{\tau_{0}}{T}-1}d\left(l\right)$. It contains exactly one period of the cyclic symbol stream and it shifts cyclically with time and, thus $p_{m}$ is constant. Accordingly we can simplify the ACF feature for periodical signals to

$$\left|\rho\left[x\right]\left(\tau_{0}\right)\right|=\left|\sum_{k=0}^{N-\tau_{0}}e^{\pi\eta i\left(-\sum_{l=j_{k}}^{j_{k}+\frac{\tau_{0}}{T}-1}d\left(l\right)\right)}\right|=\left(N-\tau_{0}+1\right).$$

The complex value of the ACF of a truely periodic signal is independent of modulation index and symbol stream.

In the second step we assume a partly-periodic signal x of ratio $p_{\%}$ with lag $\tau_{0}$. The ACF of a uniformly distributed series of samples converges to zero at all non-zero lags. Thus, according to our model the autocorrelation of the data part of the frames all but vanishes. To simplify the analysis, we introduce an indicator function b; indicating wether a sample is part of the periodic share of a signal:
$$\begin{array}{ccc}
b\left(k\right) & := & \begin{cases}
1 & k\in\text{header-part of x}\\
0 & \text{else}
\end{cases}\end{array}$$

We get

\begin{array}{cc}
\rho\left[x\right]\left(\tau_{0}\right)\\
= & \sum_{k=0}^{N-\tau_{0}}b\left(k\right)e^{i\left[\pi\eta\left(\left(\frac{k}{T}-j_{k}-\frac{k+\tau_{0}}{T}+j_{k}+\frac{\tau_{0}}{T}\right)d\left(j_{k}\right)\right)\right]}e^{i\left[\pi\eta\left(\sum_{l=0}^{j_{k}-1}d\left(l\right)-\sum_{l=0}^{j_{k}+\frac{\tau_{0}}{T}-1}d\left(l\right)\right)\right]}\\
= & \sum_{k=0}^{N-\tau_{0}}b\left(k\right)e^{i\left[\pi\eta\left(\left(-\frac{\tau_{0}}{T}+\frac{\tau_{0}}{T}\right)d\left(j_{k}\right)\right)\right]}e^{\pi\eta i\left(-\sum_{l=j_{k}}^{j_{k}+\frac{\tau_{0}}{T}-1}d\left(l\right)\right)}\\
= & \sum_{k=0}^{N-\tau_{0}}b\left(k\right)e^{\pi\eta i\left(-\sum_{l=j_{k}}^{j_{k}+\frac{\tau_{0}}{T}-1}d\left(l\right)\right)}.
\end{array}

The term $p_{m}$ (see [eq:ACF_FSK_PhaseMemory]) is not constant like it was above. As the frames are not truely periodic, the data-part of the frames introduce symbols randomly choosen from a uniform distribution. We split the results by modulation index:

If the modulation index $\eta$ is integer, $\eta p_{m}\in\mathbb{Z}$ and we get $e^{\pi\eta i\left(-\sum_{l=j_{k}}^{j_{k}+\frac{\tau_{0}}{T}-1}d\left(l\right)\right)}=1.$

Thus,
$$
\left|\rho\left[x\right]\left(\tau_{0}\right)\right|=p_{\%}\left(N-\tau_{0}+1\right).$$

If $\eta$ is not integer,

$$\left|e^{\pi\eta i\left(-\sum_{l=j_{k}}^{j_{k}+\frac{\tau_{0}}{T}-1}d\left(l\right)\right)}\right|=1$$

but the argument of the summands in ${arg}\left(e^{\pi\eta i\left(-\sum_{l=j_{k}}^{j_{k}+\frac{\tau_{0}}{T}-1}d\left(l\right)\right)}\right)$

is not constant, as it depends on the random data frame, and thus  consists of a sum of complex exponentials of differing phase. The sum decreases compared to the magnitude of the sum of an integer modulation index signal. I conclude that a random walk is in order to determine the effects of the cancellation.


Filtering noise from FFT, where noise is known to be shifted delta function I have an input sequence $x=\{x_1, x_2 , ... x_n\}$ of reals, where $n=2^m$ for some $m$. I wish to calculate FFT of $x$.

$X=FFT(x)$

However, before  I calculate the FFT, the signal $x$ gets corrupted with noise $\eta$, so $\hat{x}=x+\eta$, and calculated FFT, $\hat{X}$ is FFT of $\hat{x}$, rather than $x$.

$\hat{X}=FFT(\hat{x})$

I wish to  recover actual or approximate FFT $X$ from $\hat{X}$ under following assumption:

$x$ is band limited signal; FFT spectrum dies down quickly.
The noise $\eta$ is a shifted delta function, with unknown shift and magnitude, i.e. $\eta(n)= R\delta (n-l)$, where $R,\ l$ are real and integer respectively.

How will I do it in a compute efficient way?

High-pass filter the input with a zero-phase filter to eliminate the signal of interest, leaving only the high-passed noise. If the signal is periodic, you can do this in the frequency domain and convert back to time domain by IFFT. Find the time of the peak and deduce the magnitude of the unfiltered delta function from the magnitude of the filtered delta function (the peak value or root mean square). Knowing the delay and magnitude, subtract the noise from the original input signal. Do FFT.


audio reprocessing for machine learning I have a bunch of phone interviews where there is large sound volume difference between two voices in some of the interview audios, because they were recorded only on one side of the phone. May I ask what is the preferred procedure for reprocessing these audios so they are fit for downstream machine learning feature extraction and modelling tasks while the characteristics of the audio will not be sacrificed?

After some reading on audio compression, normalisation, and levelling, e.g. here, I have to say I am still not sure what to do next. I have seen simple minmax scaling being suggested for normalising the waveform amplitude. I have also been suggested to normalise the audio to a target amplitude (e.g. normalize in pydub.effects or ffmpeg loudnorm/dynaudnorm), which is loudness normalisation if I am not mistaken. Should I use a compressor or automatic gain controllor? Descript has a auto-leveling feature that seems to be useful, though I can't find any technical info on it.

My current plan, which is the simplest one, is to segment the audio streams to speaker turns, then normalise each speaker turn to (-1, 1) by minmax scaling. How does this sound? Btw any suggestions on declipping? Some of the audios have been clipped. Thanks!

The decision on which method to use to scale the input is very much determined by the objective and therefore what follows the scaling.

A simple linear scaling (whether peak, minmax or other) propagates to the rest of the processing chain as a multiplication. But anything that affects the dynamics of the signal (how quickly it rises and falls) ultimately shapes its content and might affect the classification results. I am assuming here, that this classification is based on some sort of spectral features (?)

The point is to look at the feature space and make sure that it is more or less spread uniformly across its bounds and representative of each class. Any improvement brought about by scaling of the original dataset would be assessed on this basis.

If your original material is something like a talk show, or other occassion where there is clean speech and no overlap of speakers then plain linear scaling would be fine or not even required.

But if you have a heated up discussion or other "quirky" content where people might be coming near the microphone or there are shouts or unpredictable background noise or people talking and overlapping, then that is much more difficult to resolve. A compressor could help up to a point but you may find that you have to provide more training samples to the classifier anyway, so that it learns that "Calm speaker A" and "Shouting speaker D" is the same person (i.e. scaling would not help much there).

  ...any suggestions on declipping?

In reality, if the audio clipped the information that was to be recorded at that point is entirely lost. And it should be treated as entirely lost.

But, if you have a huge amount of data and not quite but somewhat clipped audio. Maybe it was clipping but not all the time, maybe it was saturating and popping but not entirely hitting the rail and staying there for 5 seconds (for example),....then it might be possible to treat this as a regression type of problem and extrapolate what the clipped part was supposed to look like. This however is easier said than done because there are a lot of parameters that affect the performance of the declipping. For more information see here.

Hope this helps.


Identify rectangles of different colors How identify the presence or absence of regular stripes of different colors, but ranging from very very very very light pink to black inside of a scanned image (bitmap 200x200dpi 24-bit).
Carry a few examples.

Example 1

Example 2 (the lines are in all the columns except 7 in the second row of the last column)

For now try to identify (using python language) whether or not there is at least 5-10 pixels for the presence of different color from white to each strip, however, does not always work because the scanned image is not of high quality and the strip changes color very similar to color that surrounds it.

Thanks.

My suggestion would be in following step.

If we know that outer shapes (rectangles) are at regular interval - first we should do some sort of segmentation and/or localization.
Remove portions of background which are of no use and each block is identified and processed separately.
Now if you are expecting lines only horizontally, you can traverse vertically and create scan lines from top-to-bottom. If the region has no dark separations - the scan lines will appear straight albeit with noise. If it crosses the dark line the intensity in the scan line will drop. This will tell the position of dark strip in the block.
A set of 4 - 5 scan lines (or more) will provide a good co-relating answer.


How to add a fake thermal effect to images I am working on a project to convert a normal RGB image to thermal images. I applied color maps using Matplotlib library and got an array which then I converted back to the RGBA image using Image.fromarray() available in the PIL library.
Expected Output:

Actual Output

My code is given below:
from PIL import Image, ImageOps
import matplotlib as mpl
import matplotlib.cm as mtplitcm
import PIL
import numpy as np
from matplotlib import cm

input_image = Image.open("/Users/adwaitnaik/Desktop/github/blendpy/test_data_for_yolo/lion.png")
gray = ImageOps.grayscale(input_image)

gray.show()

#initialize the colormap
colormap = mpl.cm.jet

#normalize
cNorm = mpl.colors.Normalize(vmin=0, vmax=255)

#initialize the mapping
scalarMap = mtplitcm.ScalarMappable(norm=cNorm, cmap=colormap)

#create image array
image_array = scalarMap.to_rgba(gray)

#convert to numpy array
image_np = np.asarray(image_array)

image_np

#image_array
image_colored = Image.fromarray(image_np,'RGB')

image_colored.show()

This is the image I used as Input

Can anyone suggest a workaround?
Thanks

After experimenting with a couple of libraries, my code finally worked with OpenCV. Code and Output can be seen below.
Code
#this code uses color maps which works as thermal camera filter on the input image

import sys
import numpy as np
import cv2
import matplotlib as mpl
import matplotlib.cm as mtpltcm

def main(argv):

    while (True):
        # read the file
        input_img = cv2.imread('path_to_the_file')

        # convert image into grayscale
        gray = cv2.cvtColor(input_img, cv2.COLOR_BGR2GRAY)

        #initialize the colormap
        colormap = mpl.cm.jet

        #add a normalization
        cNorm = mpl.colors.Normalize(vmin=0, vmax=255)

        #init the mapping
        scalarMap = mtpltcm.ScalarMappable(norm=cNorm, cmap=colormap)

        colors = scalarMap.to_rgba(gray)

        #to display the image
        cv2.imshow("_animal_name",colors)

        if cv2.waitKey(1) & 0xFF == ord('q'):
            break

    cv2.destroyAllWindows()

if __name__ == '__main__':
    sys.exit(main(sys.argv))

Expected Output


Smoothing complex data by convolution I need to smooth noisy complex data with a Gaussian filter.

Right now, I apply the filter to real and imaginary part of the data separately, which needs two convolutions. The intended results are the equally smoothed components of the array.

In order to reduce computational load, would it be possible to filter the complex array as a whole and take real and imaginary part after smoothing?

If you filter complex data $x[n]$ with a filter with a real-valued impulse response $h[n]$, your complex-valued output is given by the convolution sum

$$y[n]=\sum_kx[k]h[n-k]\tag{1}$$

The real and imaginary parts are

$$y_R[n]=\sum_kx_R[k]h[n-k]\tag{2}$$

and

$$y_I[n]=\sum_kx_I[k]h[n-k]\tag{3}$$

So whatever you do, in order to compute the complex-valued output $y[n]$ you need to perform two real-valued convolutions / filtering operations.


Why is histogram equalization called as "equalization"? When output histogram is not same as reference histogram, why is the method called histogram equalization?

For example below are two figures from http://en.wikipedia.org/wiki/Histogram_equalization. The output histogram (2nd figure) looks nothing like a discrete uniform distribution. Why is it called "equalization" then?

To answer your question histogram equalization is called like this because its function is to produce an equalized histogram (that is an uniform probability density function).

There are different algorithms that may approach this function, and obviously there is a problem in the example that is shown:

In fact, the algorithm used there will always have trouble producing a flat histogram when mapping the function to a discrete set of 256 (except when the input histogram is already flat!)
As a consequence, the cumulative histogram shows discrete jumps which produce this non-regular output histogram.

Why is this simple algorithm still used? Because,

the output image will often good enough for the initial purpose,
that is to maximize the range of intensities used in the image - though the histogram may look ugly...
it is simple!

There are simple alternatives. The simplest is to use a linear interpolation of the cumulative function.

Theory: In general terms, the method may be easily derived if we know the probability distribution function $dP_{i}$ of variable $a_i$ representing the luminance image at pixel $i$. Indeed,  by choosing the non-linearity as the cumulative distribution function transforming any observed variable $\bar{a}_i$ into:
$$
z_i(\bar{a}_i)= P_{i}(a_i \leq \bar{a}_i) = \int_{-\infty}^{\bar{a}_i} dP_{i}(a_i)%
$$
This is equivalent to the change of variables which transforms the luminance image to a variable with uniform probability distribution function in $[0, 1]^M$.
Evaluation: In physiology (retina for instance), this distribution is assumed to be sampled for each "pixel" (photoreceptor in mammals) over time (for instance, something like "the probability distribution over last minute"). In image processing, when handling one single image, you can assume stationarity over space, that is compute the histogram over all pixel values.
Algorithm: now that we have all in place, some example algorithm in python:

Sample code:

import numpy as np
from scipy.misc import lena
I = 1. * lena()
I = (I - I.min())/(I.max()-I.min())
import matplotlib.pyplot as plt
# the cumulative histogram
n_quant = 40
(hist, bins) = np.histogram(I.ravel(), bins=np.linspace(0., 1., n_quant+1), normed=True)
P = np.cumsum(hist/np.sum(hist))
# equalization
I_z = np.interp(I, np.linspace(0., 1., n_quant), P)
(hist_z, bins) = np.histogram(I_z.ravel(), bins=np.linspace(0., 1., n_quant+1), normed=True)
plt.bar(np.linspace(0., 1., n_quant), hist, color='r', alpha=.3, width=1./n_quant)
plt.bar(np.linspace(0., 1., n_quant), hist_z, color='b', alpha=.3, width=1./n_quant)
plt.title('Histogram of Lena')
plt.xlabel('luminance')
plt.ylabel('probability')
plt.xlim((0,1.))
plt.savefig('equalization.png')

The output (red is raw histogram, blue is equalized). As you can see, there are variations, but it looks better :-). Notably, you have less "jumps" around most popular values. The parameter n_quant is important to tune smoothness versus precision.


segmenting region of interest (ROI) from background I am dealing with images from biological slides and want to segment out the ROI (cells) from the background. I have tried k-means clustering, otsu thresholding, adaptive thresholding, edge detection etc. Running out of ideas now.
Which other segmenting techniques can be tried for such images? Keeping in mind that the color of cells and background can change
click here for image

Some very cool although a bit more time consuming options:

Active contours will try to find the best "edge" for the cells while maintaining smoothness.
and is implemented in matlab2013a ( active contour)
You might try using mser which is an approach which detects regions which are similiar and in color and stable over the image. which means they stay about the same size over a wide array of thresholds.
I used it for a similar problem with detecting ultrasound tumors.
If you have a good database of images ( good = covering most cases )  I would try some machine learning algorithms. for example, you could learn on the mean and variance of the objects in the image.
If you have trouble separating the cells after segmentation ( i.e. you have two connected cells ) I would try looking at the watershed function.


Analytical Solution for the Convolution of Signal with a Box Filter I have an exercise in which I am trying to filter an input signal $y(x) = \sin(x)$. Ideally, I would like to apply a box filter to this signal.

Previously, I successfully convolved the input signal $y(x)$ with a decaying response $h(x) = e^{-x}$.

I did so by the following the definition of convolution (e.g., integrating $\int_0^t\sin(x')e^{-(x-x')}\mathrm{d}x'$ and computing a damped sinusoidal signal.

My box filter is given by $\frac{1}{\Delta}$ for$|x-\xi| \leq \frac{\Delta}{2}$ and 0 elsewhere, where $\Delta$ is the filter width. I understand that a box filter is a local average, and I can implement this numerically, but I do not understand how to analytically integrate this as I did with the damped exponential 'filter'.

I tried to take the Fourier transform of $y(x)$ and $h(x)$ and multiplying them in Fourier space, but I could not figure out how to do so.

Thanks for any help.

DO NOT PASS GO; DO NOT Collect $200; DO NOT Take Fourier transforms, or worse yet, FFTs

You do the convolution exactly the way you would do any other convolution: start with the basic convolution integral (not what you wrote) and apply the properties of the signals that you are using to come up with an easier calculation.  Begin with
$$\int_{-\infty}^\infty y(x^\prime)h(x^\prime-x)\mathrm dx^\prime\tag{1}$$
and use the result that $h(\cdot)$ is nonzero only when its argument ($x^\prime-x$ in this instance) lies in the interval $\left[-\frac{\Delta}2, +\frac{\Delta}2\right]$. So, the integrand of the convolution integral is $0$ whenever $x^\prime$ is such that
$$x^\prime-x > \frac{\Delta}2 \implies  x^\prime > x + \frac{\Delta}2$$
or that
$$x^\prime-x < -\frac{\Delta}2 \implies  x^\prime < x - \frac{\Delta}2.$$
This allows us to simplify the convolution integral $(1)$ into
$$\int_{x - \frac{\Delta}2}^{x + \frac{\Delta}2} y(x^\prime)\frac{1}{\Delta}\mathrm dx^\prime\tag{2}$$
where I have substituted the nonzero value of $h$ in the integral to save a step.  Can you take it from here?


Calculating signal/noise ratio and the noise level If the received signal strength for a 5.2 dBi dipole was measured at -46 dBm and the noise level was -96 dBm, what would the signal/noise ratio be in dB, the received power in µW, and the noise level in nW?

The decibels (dB) provides a logarithmic scale for the ratio of two power values. Mathematically it is defined as:
$$
  R = 10 \log_{10} \left(\frac{P}{P_0}\right)
$$
where $P$ and $P_0$ are power values in linear scale (such as Watts) and $R$ is the ratio expressed in logarithmic scale. Often a fixed reference power level (such as 1 Watt) is chosen for $P_0$. When that reference power level $P_0$ is 1 milliwatt (mW), the resulting unit for the $R$ quantity is called a dBm.

Given a known reference level, you can conversely obtain the power in linear scale from the power ratio value expressed in logarithmic scale by using the inverse formula:
$$
  P = 10^{R/10}\cdot  P_0
$$

So for a measurement of -46dBm (remember that the dBm units means that the reference power level is 1mW), the resulting power is:
$$
\begin{align}
  P_s &= 10^{-46/10} \cdot 1 \mbox{mW} \\
      &\approx 2.5 \times 10^{-5} \cdot 1 \mbox{mW} \\
      &\approx 2.5 \times 10^{-8} \, \mbox{W} \\
      &\approx 2.5 \times 10^{-2} \, \mbox{µW}
\end{align}
$$
Similarly for a measurement of -96dBm, the resulting power is:
$$
\begin{align}
  P_n &= 10^{-96/10} \cdot 1 \mbox{mW} \\
      &\approx 2.5 \times 10^{-10} \cdot 1 \mbox{mW} \\
      &\approx 2.5 \times 10^{-13} \, \mbox{W} \\
      &\approx 2.5 \times 10^{-4}  \, \mbox{nW}
\end{align}
$$

Now for the signal-to-noise ratio (SNR), it is simply defined as the ratio of the signal power to the noise power. So in linear scale, that would be:
$$
\begin{align}
  \mbox{SNR} &= \frac{P_s}{P_n} \\
             &\approx \frac{2.5 \times 10^{-8}}{2.5 \times 10^{-13}} \\
             &\approx 10^{5}
\end{align}
$$
which when expressed in decibels would give us:
$$
\begin{align}
  \mbox{SNR}_{\mbox{dB}} &= 10\log_{10} \left(\frac{P_s}{P_n}\right) \, \mbox{dB} \\
             &= 10\log_{10} \left(10^{5}\right) \, \mbox{dB} \\
             &= 50\mbox{dB}
\end{align}
$$
We could of course also have taken a bit of a shortcut to compute that SNR directly from the power measurements in dBm by using logarithmic identities to get:
$$
\begin{align}
  \mbox{SNR}_{\mbox{dB}} &= 10\log_{10} \left(\frac{P_s}{P_n}\right) \, \mbox{dB} \\
             &= 10\log_{10} \left(\frac{\frac{P_s}{P_0}}{\frac{P_n}{P_0}}\right) \, \mbox{dB} \\
             &= 10\log_{10} \left(\frac{P_s}{P_0}\right) - 10\log_{10}\left(\frac{P_n}{P_0}\right) \, \mbox{dB} \\
             &= \left(-46\,\mbox{dBm}\right) - \left(-96\,\mbox{dBm}\right) \\
             &= 50\mbox{dB}
\end{align}
$$


How to account for multiple signals in the input to a Kalman filter I am interested in using a Kalman filter to track an oscillating signal due to a physical oscillator in a signal which contains the signal from multiple such oscillators oscillating at different frequencies. The signal in this case is simply a voltage value changing with time.

I have thus-far been using sinusoidal Kalman filters for this, however they do not perform optimally as the oscillator is not a perfect sinusoid, it experiences amplitude changes, phase jumps and changes in frequency slowly over time, the first 2 it can correct for, but the last causes the Kalman filter to not be centred on the correct frequency anymore and perform poorly.

However it was suggested to me, that, as long as I use a high sample rate (~25 samples per oscillation) a simple kinematic filter (tracking position, velocity and acceleration) will track the signal, and then there is no problem if the frequency shifts as it is just tracking the position of the oscillator with no forced frequency assumed.

For this kinematic filter, would it require me to filter out the signal of the oscillator to be tracked from the other oscillators with a bandpass filter? Alternatively would tracking the motion of all the oscillators simultaneously work? Without some sort of technique I see no way that such a filter would know what the signal it was required to track was.

The kinematic filter (which is mathematically a polynomial Kalman filter) has the ability to track any signal as long as the sampling rate is high enough, because when it's so, the complex motion of the target becomes linearized and fits into the quadratic kinematic model being employed by the Kalman filter dynamic model but the disadvantage is the processing cost due to high sampling rate.

Also the polynomial (kinematic) Kalman filter will not distinguish between two sinusoids but instead track their sum ignoring the noise only, unlike a specifically designed sinusoidal filter which tracks the chosen sinusoid and ignores all the rest as disturbance. However the model needs to be exact and known in advance as its fundamental disadvantage. The figure below shows a simulation example of a sinusoidal Kalman filter with two sinusoidal signals added together plus uncorrelated zero mean white noise as the measurement signal. It can be seen that the sinusoidal Kalman filter tracks the chosen sine wave (which is used in deriving the system dynamics matrix) and ignores all the rest as disturbance including the added white noise and the other sinewave which was at an harmonic multiple of the tracked sinusoidal frequency.

The extended kalman filter is a compromise, it has some additional computational burden, specifically the per sample update of the linearized derivatives, but a compromise can be found by lowering the sampling rate compared to a kinematic filter's requirement.


Image zoom out using Gaussian blur and downsampling I am trying to implement a simple zoom out algorithm using two steps:

The image is blurred using gaussian convolution of standard deviation $0.8\sqrt{{\tt zoomfactor}^2-1}$ as we saw in class.
The image is then downsampled using a simple downsampling of zoom_factor step
The following code is supposed to do the trick.

In class, we saw that zoom out using only downsampling should give worst results than in the case of zoom out combined with gaussian blur. However for my case, this was the opposite which is really weird.

 close all;
 clear all;
 I = double(imread('put an image path here'))/255;

 zoom_factor=2;
 %I=rgb2gray(I);
 sigma=0.8*sqrt(zoom_factor*zoom_factor-1);
 [M,N,s]=size(I);
 f1=-fix(M/2):ceil(M/2)-1;
 f2=-fix(N/2):ceil(N/2)-1;

 [f2,f1]=meshgrid(f2,f1);
 Y=exp(-2*pi*pi*sigma*sigma*((f1/M).^2+(f2/N).^2));
 If=fft2(I);
 If=fftshift(If);
 for i=1:3
   If(:,:,i)=If(:,:,i).*Y;
 end
 newI=ifftshift(If);
 newI=ifft2(newI);

 newI=newI(1:zoom_factor:end,1:zoom_factor:end,:)

 imshow(newI)

What I understood from your post is that you don't like the excess blur that happens at the downsampled image.

The theory of signal processing states that, when you downsample a digital signal by a factor of $M$ there's the potential of aliasing to happen if the signal is not bandlimited to $|w| < \pi/M$;

The spectral effect of the downsampling is such that the non-zero portions of the original signal spectrum above $ |w| > \pi/M$ can overlapp into fake frequency positions in the range $-\pi< w<\pi$ so that they will create visual artefacts that does not belong to the original signal.

In order to avoid such an overlapp (to avoid aliasing) the formal approach is to apply a lowpass filter with a cutoff frequency of $w_c = \pi/M$ and a gain of 1 at the passband.

The effect of this filter is to blur the image to the extend that when downsampled there won't be any significand overlap from the high frequency portions of the image. On the other hand this is eventually a blurring operation and it will throw away signal content (details of the image). So there will be a loss of quality (in terms of high frequency details) in your downsampled image.

However the problem in your code is that there is too much blurring, due to the Gaussian filter radius set to an inappropriate value. In theory what you need is an ideal filter with sharp transition at the cutoff frequency $w_c = \pi/M$ but the filters desiged in your code produces not only mild transition region but also a lower cutoff frequency than necessary.

So you can modify your lowpass filter radius untill you get an acceptable downsampled image.

Below is the modified code that shows spectrum of the Gaussian filter, you can use this to deduce if the filter cutoff is too much.

close all; clear all; clc;
I = double(imread('Cameraman.tif'))/255; % Black and Whit eimage in [0-255] range
figure,imshow(I);title('Original Image');

zoom_factor=2;

sigma=0.8*sqrt(zoom_factor-1);   % I modified the radius here ???
[M,N,s]=size(I);
f1=-fix(M/2):ceil(M/2)-1;
f2=-fix(N/2):ceil(N/2)-1;
[f2,f1]=meshgrid(f2,f1);
Y=exp(-2*pi*pi*sigma*sigma*((f1/M).^2+(f2/N).^2));  % The BLUR filter
figure,imshow(abs(Y));title('Spectrum of the filter')

If=fft2(I);
If=fftshift(If);
for i=1:1
  If(:,:,i)=If(:,:,i).*Y;
end
newI=ifftshift(If);
newI=real(ifft2(newI));
figure,imshow(newI);title('Gaussian Filtered Image')

newI=newI(1:zoom_factor:end,1:zoom_factor:end,:);
figure,imshow(newI);title('filtered and downsampled')

figure, imshow(newI);title('filtered and downsampled image')

figure, imshow(I(1:zoom_factor:end,1:zoom_factor:end));
title('downsampled without filtering')

Another solution is to use a highpass filter after the downsampling operation to improve the image by emphasizing the high frequency detail a bit.


Area under waveform in t-domain including 2 delta functions question:

my attempt:

$(a)$ $\displaystyle\int_{-\infty}^{\infty} f(t) dt=2 u (t)+r(t)-r(t-T)+Ku(t-T)$

then we draw it's diagram and found it was unbounded so,we doubt how area can be zero for any value of $K$ so, we need help here finding value of K (which we think  will come negative )

$(b)$ $\displaystyle\int_{0_{+}}^{\infty} f(t) dt=r(t)-r(t-T)+Ku(t-T)$

same diagram as part $(a)$ , and same story as above ,except  we pull above curve down by two units .

so , we're unable to find value of $K$. any help will be appreciated

thanking you !

edit : $r(t)$ is ramp function $= t u(t)$

please someone confirm if following  devlopment is right or  wrong .becuase answer of exercise question isn't provided in book "F F kuo"

we know

$f(t)\xrightarrow{\mathcal{F}}F(j\omega); $ where, $F(j\omega) =\displaystyle\int_{-\infty}^{\infty} f(t) \ e^{-j\omega t} \ dt\implies F(0)=\text{area under time domain curve }$

So, $(a)\ F(j\omega)=2+T.\ Sa{(\omega T/2)}+Ke^{-j\omega T}\implies F(0)=2+T+K=\displaystyle\int_{-\infty}^{\infty}f(t)dt$

but area is given zero therefore , $K=-(2+T)$ .

and it is resonable also becuase strength of impulse at $t=T$ must be equal to the negative summation of  area of regions preceding it.

similarly,

for part $(b)\ $   $K=-T$ because time is starting from $t=0_{+}$ that's why excluding impulse of strength  $2$ at $t=0$

note :

$Sa (\omega T/2)=\dfrac{\sin\left(\dfrac{\omega T}{2}\right)}{\omega T/2}$


MATLAB and Radon Transform-cannot obtain a sinogram I'm quite new to MATLAB and I'm trying to obtain an implementation of the computed tomography without using any built-in functions such as radon() and iradon(). I guess I understood what Radon Transform does: it simply adds up projections taken at some angles (in matrix columns in MATLAB) and produces a sinogram as a result.  But I am at a loss at this point. My attempt never produces anything that is similar to a sinogram. What am I doing wrong? I cannot get ceil(128*sqrt2) as the maximum value for the detector size as well (of one unit each on the projection axis).

Side note: Well, I admit that the code is lousy (my image, too), but I tried my best in the last week to decipher MATLAB. Thank you very much in advance for any kind of help and suggestions.

A=zeros(128,128);
for k=1:128
    for n=1:128
        if k>70&&k<100&&n>70&&n<100
            A(k,n)=1;
            else if k>20&&k<35&&n>30&&n<55
            A(k,n)=0.7;
                else if k*k+n*n>100*40&&k*k+n*n<100*50
                    A(k,n)=0.8;
                    else if (k-50)^2+(n-50)^2<=400
                        A(k,n)=0.5;
                        else if (k-100)^2+(n-20)^2>=100&&(k-100)^2+(n-20)^2<=400
                                A(k,n)=0.8;
                            end
                    end
                end
            end
        end
    end
end
imshow(A)
M=size(A);
[ilength, iwidth]=size(A);
idiag=sqrt(128^2+128^2);
lengthpad=ceil(idiag-ilength)+2;
widthpad=ceil(idiag-iwidth)+2;
pad=zeros(ilength+lengthpad,iwidth+widthpad);
pad(ceil(lengthpad/2):(ceil(lengthpad/2)+ilength-1),ceil(widthpad/2):(ceil(widthpad/2)+iwidth-1)) = A;
figure, imshow(pad)
N=100;

% angles through which we will perform the projection
for j=0:N-1
            thetaj=j*(pi/N);
            costhetaj=cos(thetaj);
            sinthetaj=sin(thetaj);

    %finding the midpoints of the pixels on the image
    for iy=1:size(pad)
        for i=1:size(pad)

            y_img=iy-pad/2-0.5;
            x_img=i-pad/2-0.5;
            img=sum(costhetaj*x_img*pad(i,iy)+sinthetaj*y_img*pad(i,iy),1);
            last_img(i,iy)=sum(img);

        end
        end
end
    figure, imshow(last_img)
    figure, plot(last_img)

Edit: I added zero padding to the code, but I still don't know why I get half of the image as output.

There are two points that stand out from the provided code:

The determination of detector_size

You seem to be trying to track the "widest line of sums" (?). But we know that this is the diagonal and it occurs at 45 degrees in square images. Therefore, we can "pre-allocate" our detector to have that maximum size anyway.

The rotation of the image

The code seems to be along the right track, taking the "image grid coordinates" (i,iy) and rotating them to (x_img,y_img) but:

There is no centre of rotation in the calculation (?). This usually is taken to be at the centre of the image.
While some rotation seems to be taking place, the data is still read from the original locations of A (involving some scaling too) (?).

There are two approaches that you can take here:

Rotate the image
Rotate the "detector" (but do it properly)

The first approach is dead simple but a bit wasteful in terms of memory and computation. It goes a little something like this:

Given a $M\times N$ image ($I$), construct a $D \times D$ image ($A$) where $D = \sqrt{M^2 + N^2}$
Now, centre $I$ into $A$. This basically translates $\left(0,0\right)$ of $I$ to $\left(\frac{D-M}{2},\frac{D-N}{2}\right)$ in $A$.
For each  angle $\theta \in \Theta$, rotate the image (using imrotate) to produce $A_r$. Now, note here that depending on your frame of reference, you might have to transform $\theta$. That is because, if you sum the columns of the rotated image, then your detector is positioned at $+ \frac{\pi}{2}$. But that is just a simple linear transformation to align your results with convention. Also important to note here that it is worth rotating $A$ afresh each time rather than keep rotating $A_r$ by increments, which would progressively distort the image due to successive applications of interpolation.
Once you have $A_r$, produce the sum of the "detector" by adding up its columns.

And that is it.

It is dead simple because you don't have to handle the rotation of the image which, by itself, is not difficult but it involves interpolation which can get messy. The rotated image grid will never land on integer pixel locations and so you are going to have to apply some form of interpolation that goes beyond nearest neighbour. For linear interpolation (for example), this would involve determining the patch of the 4 nearest pixels to the target location, then obtaining the distances to those pixels and finally creating the weighted sum.

It is "wasteful" in CPU because $A$ is bigger than $I$ and all pixels are rotated even if their contribution to a sum is in fact zero.

It is "wasteful" in memory because $A$ is bigger than $I$ and for large differences of $M,N$, large square images are produced.

The immediate improvement you can apply to this is to rotate the "scanning grid" (as you are trying to do now) but when you obtain a slice off of the image to centre it in the "detector slice" which would have to be $D$ pixels wide anyway. In this case of course, you still have to handle interpolation. But in this approach, any coordinates landing outside of the original image are simply skipped.

Hope this helps

EDIT:

The rotation is still not performed correctly:

% Load the image
Q = imread("testImage.jpg");
% Resize the image, in case it is too big
Q = imresize(Q,0.05);
Q_size = size(Q);
% Find its diagonal
D = ceil(sqrt(sum(Q_size.^2)));
% Create a square image of size equal to the diagonal
I = zeros(D,D);
I_size = size(I);
% Center the image in the square
offset = ceil((D-Q_size)./2);
% Copy the original image across
I(offset(1):offset(1)+Q_size(1)-1,offset(2):offset(2)+Q_size(2)-1) = Q;

I contains the padded image.

Using imrotate:

% Theta angles to obtain sums at
theta = 0:11.25:180;
% Initialise the output. One row for each angle, one column for every summing "sensor".
R1 = zeros(length(theta), D);
for k=1:length(theta)
    R1(k,:)=sum(imrotate(I,theta(k),'nearest','crop'));
end;

figure;imshow(R1./max(max(R1)))

The above finishes very quickly. Notice nearest as interpolation in imrotate and crop so that imrotate does not add yet more padding.

Not using imrotate:

% Rotating the detector
theta_rad = (theta./180).*pi; %This simply converts theta (from above) in radians
% Initialise the output. One row for each angle, one column for every summing "sensor".
R2 = zeros(length(theta_rad), D);
% THE CENTRE OF THE IMAGE
I_centre = round(I_size./2);
% FOR EVERY ANGLE
for k=1:length(theta_rad)
    % m,n scan the image left to right top to bottom...BUT...
    for m=0:(I_size(2)-1)
        for n=0:(I_size(1)-1)
            % THEY GET REMAPPED TO U,V
            % Remap m,n to their rotated versions.
            % I do the rotation in polar coordinates here
            rho = sqrt(sum((I_centre-[n,m]).^2));
            phi = atan2(n - I_centre(1), m - I_centre(2));
            % This point now maps to
            u = round(I_centre(2) + rho .* cos(phi-theta_rad(k)));
            v = round(I_centre(1) - rho .* sin(phi-theta_rad(k)));
            % Check that the point is inside the image
            if (u>0 && u<I_size(1) && v>0 && v<I_size(2))
                % NOTICE THAT I IS READ FROM THE REMAPPED LOCATIONS
                R2(k,m+1) = R2(k,m+1) + I(v,u);
            end;
        end;
    end;
end;

figure;imshow(R2./max(max(R2)))


4G LTE Sampling Frequency Offset Could someone please explain how the 4G LTE OFDM measures and corrects the sampling frequency offset?
I can see a number of algorithms proposed in the literature, but I am trying narrow it down to the actually implemented ones.
Thanks in advance.

I am not a specialist of the 4G LTE, but worked on other OFDM waveform similar to it.
OFDM signal is made of symbols.
Each symbol in frequency domain is made of a number of subcarriers.
Some of subcarriers are mapped with data and other with pilots.
Pilots have known defined values (usually scrambled BPSK symbols).
The sampling frequency offset results in phase shift of the same pilot from one OFDM symbol to the next one (or next ones).
After computing the phase shift, it can be used as a timing error of a digital phase locked loop acting on the resampler of the demodulation system.
You could try to analyze the behavior of pilots with sampling offset if you have a MAtlab model.


How to find SPL of a WAV file given a calibration tone in WAV format My ultimate aim is to find the Sound Pressure level of a digital WAV audio file and represent this (say, in MATLAB) as a dB SPL re 20 uPa on y-axis and frequency on x-axis.

As this link mentions:

  Calibration Signal Recording – Record a WAV file with a tone or signal
  of a known level using the same recording input gain setting that were
  used to record the WAV file.  Sound Quality Analyzer will use the
  calibration recording to generate the proper normalization term for
  the WAV file to be analyzed.

Now, let'say I have a calibration tone which I know has a 94 dB SPL at 1 KHz frequency, how can I use this information to scale the pressure values my audio file to finally get to a dbSPL.

In the above equation, I know Po = 20 uPa (threshold of hearing). So if I am able to find the scaled P values of the audio file, I can use this equation to find the dBSPL.

Edit:

So the question boils down to this:

There is a noise.wav that needs to be analysed.
There is a calibration.wav that is recorded with 94 dB SPL at 1 KHz. This is the known calibration tone from a calibrator. The settings on the calibrator say 94 dB, 1KHz. Is this 94 dB SPL re 20 uPA? How do we find the SPL of this recording to get a normalization factor for the noise.wav to be analyzed? I am assuming the actual SPL of the recording will be different than what the calibrator says.
How will you use the normalization factor obtained above to calibrate and find the SPL of the original noise.wav? If noise.wav is a 1 minute recording - do we need to take n number of samples (window the signal), calculate the SPL re 20 uPa of each window and further normalize the SPL value obtained from the term obtained above?

Here, I am assuming that the WAV files (noise and calibration tone) are not recorded with a calibrated sound level meter but rather, say, a mobile phone.  So, essentially what I am trying to get at is to replicate the functionality of calibrating a sound level meter given these two recordings.

The calibration file and the file you want to analyze should be captured on the same equipment and the same conditions. A reference will be obtained in a free space environment, no echoes, and far field. If you have a room impulse response, you can usually do a reasonable correction but the collection equipment should be nearly the same.
Most cheap SPL meters are usually A weighted.

There are also different ways conventions on how to report SPL. In Underwater Acoustics, we usually report dB ref 1 micro pascal in a 1 Hz band, so you need to correct your levels by the ENBW of the filter (fft bin)

There's a section in Kinsler and Fry's Fundamentals of Acoustics on converting levels for different frequency bands. Noise adds as 10log and signal adds as 20 log. I recommend you look at the book.


Laplace of step and integration are same? Why do we have Laplace transform of a step function and integrator is same.

\begin{align}
\mathcal L\left[u(t)\right] &= \frac 1s\\
\mathcal L  \left[  \int dt\right]  &= \frac 1s
\end{align}

Please clear my doubt on this.

This is because the impulse response of an integrator is $h(t)=u(t)$. The output which is the convolution with the impulse respoponse is
$$y(t)=\int_{-\infty}^{\infty}x(\tau)h(t-\tau)d\tau$$
and with $h(t)=u(t)$ it becomes
$$\begin{align}
y(t)&=\int_{-\infty}^{\infty}x(\tau)u(t-\tau)d\tau\\
&=\int_{-\infty}^{t}x(\tau)d\tau\tag{1}\end{align}$$
where $(1)$ is resulted from the fact that $$u(t-\tau)=\begin{cases}0& \forall \tau>t\\
1 & \text{otherwise}\end{cases}$$
The transfer function (which is the Laplace transform of impulse response) is $$H(s)=\mathcal{L}\{h(t)=u(t)\}=\frac{1}{s}$$


GNU Radio | Async FEC encoder/decoder loopback | Issue during decoding? I have implemented a loopback flowgraph intending to perform the following:

FEC encode (repetition code) an input stream in16_2.dat
Store the encoded stream as out16_2_fec_async_enc.dat
FEC decode the encoded stream and store it as out16_2_fec_async_dec.dat

Follows the implemented flowgraph:

My expectation: I was hoping that out16_2_fec_async_dec.dat would be equal to in16_2.dat.

What really happened: out16_2_fec_async_dec.dat was empty.
I can confirm that the encoded data, out16_2_fec_async_enc.dat, is apparently correct, by running:

$ xxd -b in16_2.dat
00000000: 00000000 00000000 11111111 11111111 01010101 01010101  ....UU
00000006: 10101010 10101010

$ xxd -b out16_2_fec_async_enc.dat
00000000: 00000000 00000000 00000000 00000000 00000000 00000000  ......
00000006: 00000000 00000000 11111111 11111111 11111111 11111111  ......
0000000c: 11111111 11111111 11111111 11111111 00001111 00001111  ......
00000012: 00001111 00001111 00001111 00001111 00001111 00001111  ......
00000018: 11110000 11110000 11110000 11110000 11110000 11110000  ......
0000001e: 11110000 11110000                                      ..

Also, the confirmation that out16_2_fec_async_dec.dat is empty:

$ xxd -b out16_2_fec_async_dec.dat

I wonder what may be missing in this setup?
Please let me know if additional information is needed.
Thank you in advance!

The issue is that the FEC Async Decoder expects the input to be soft decisions with one decision per item instead of packed PDUs.

Here's a working example:

Encoder deployment:

Packed: yes
Rev. unpacking: yes
Rev. packing: yes

Repack Bits:

Packet alignment: Output
Endieness: LSB

Decoder deployment:

Packed: yes
Rev. packing: yes


How to understand relationships between ellipse and second moment matrix of Harris corner detector? guys. I am really lacking in knowledge of linear algebra. I am reading slides of Harris corner. But I am really confused about one of them. I know that I can find corners by two large eigenvalues but I can not consider it as an ellipse....can anybody give me any suggestions?

Consider a 2-by-2 diagonal matrix $M$ with diagonal elements $a>0$ and $b>0$. Then the equation

$$x^TMx=const\tag{1}$$

(with $x$ a column vector with elements $u$ and $v$) becomes

$$au^2 + bv^2 = const$$

which is simply the equation of an ellipse (or a circle if $a=b$). If $M$ is not diagonal, equation (1) still describes an ellipse if $M$ is positive definite (i.e. it has positive eigenvalues), but then its major and minor axes do not coincide with the Cartesian axes.

Also have a look at this.


Discrete Fourier Transform of PMF and PDF I got this explenation:

"It is well-known that an image’s histogram is essentially
the probability mass function (pmf) of the image (only
differing by a scalar). Multiplying each component of the
pmf by a correspondingly shifted unit impulse results in the probability density function (pdf). Obviously, in the context
of discrete Fourier transform (DFT), the unit impulses can be
ignored, implying that we can treat pmf and pdf
exchangeable. Thus, the pdf can be thought as the
normalized version of a histogram."

and I am trying to fully understand it.

It is clear to me that histogram of an image are probability mass functions. How can we treat is as probability density function since it is discrete signal, not continuous?
Is the only reason that DFT of both PMF and PDF would give us same result so it doesn't really matter in this scenario?

And how can we even get PDF from PMF?
I don't get this part: "multiplying each component of the
pmf by a correspondingly shifted unit impulse results in the probability density function"
and I cannot find any good example or explenation.

A random variable (RV) doesn't need to be continuous in order to have a PDF. The PDF of a discrete RV, i.e. a RV that takes on only discrete values, is of course not a continuous function but is only different from zero at the possible discrete values of the RV.

An example would be the outcome of a fair die. It takes on the values 1 to 6 with equal probability $1/6$. It's PDF is given by

$$p(x)=\frac16\sum_{i=1}^6\delta(x-i)\tag{1}$$

where $\delta(x)$ is the Dirac delta impulse. The Dirac delta impulses are necessary in the PDF (as opposed to the PMF) because probabilities are computed from the PDF by integration. Since for any interval $[a,b]$ containing the value $0$

$$\int_{a}^b\delta(x)dx=1$$

the PDF given by (1) integrates to $1$, as desired. And the probability of the RV $X$ being inside an arbitrary interval $[a,b]$ is

$$P[a\le X\le b]=\int_a^b p(x)dx$$


Play with a Gaussian Random Set in the Frequency Domain to Obtain Desired Effect in the Time Domain Please assume we have a set of 100 random numbers obeying Gaussian PDF in time-domain. Let's index them 1-100.

Within a time-accurate simulation, I apply four operations on this dynamic set every time-step:

Use the entire set for some mathematical operation
Then, remove, say, the last 10 numbers of the set (indices:91-100)
Shift the remaining 90 numbers forward from the indices:1-90 to 10-100
Generate 10 new Gaussian random numbers and put them into the emptied indices:1-10. And the entire set should still obey Gaussian PDF.

I want to obtain these time-domain changes via performing some operations in the frequency domain.

However, I do not know how I should play with the frequency-domain counterpart of the set, so that I obtain the same (or similar) results as if performing time-domain operations above.

The obscure form of the question:

One can generate random samples of Gaussian distribution directly in the frequency domain as explained here, and exemplified with a Python snippet here.

Let assume we generate a one-dimensional set consisting of uniformly-spaced 100 Gaussian random samples, $X(f)$.

I want to manipulate the set as such its time-domain counterpart, i.e. $x(t)$, is:

Index-shifted by N positions in a circular manner,
Then, (only) its first N elements being replaced by new N random samples.

For the first item, this gives the answer. For the second, however, I am perplexed.

Could you please tell us how to insert this new N random samples directly into the existing set in the frequency domain, so that we can avoid possibly(?) redundant FFTs otherwise?

You're asking to do a localized operation in time using the Frequency Domain.

It's going to be not elegant, really.

Here what you can do:

Define the input signal in frequency domain as $ X \left[ k \right] $.
Multiply by Linear Phase in Frequency Domain to apply circular shift in time domain.
This means samples 1-90 will become 11-100 and 91-100 will become 1-10. Lets' call this multiplied signal $ \hat{X} \left[ k \right] $.
Pre calculate the DFT of a signal which is 1 for the first 10 samples and then zero for samples 11-100. Lets call it $ H \left[ k \right] $. This is basically a windows on the first 10 samples. In order to extract the 10 first samples in time domain we would multiply by it. In the frequency domain it means we will apply Convolution. The result of the convolution will be a frequency domain signal which is equivalent of the DFT of only the 10 first samples padded with 90 zeros. Keep this result on side. Let's Call it $ Y \left[ k \right] = \hat{X} \left[ k \right] \ast H \left[ k \right] $.
Generate noise on frequency domain which is with the properties needed of generate what's needed on time domain. Convolve the generated noise with $ H \left[ k \right] $ which will again extract a signal which is equivalent to only 10 first samples in time domain. From this signal, lets call it $ Z \left[ k \right] $ subtract the signal $ Y \left[ k \right] $.
Add the result to $ \hat{X} \left[ k \right] $. Now the signal $ \hat{X} \left[ k \right] + Z \left[ k \right] - Y \left[ k \right] $ is what you're after.

Pay attention that doing this in Fourier domain might be much less efficient than doing it in time domain.
As you need to apply convolution twice instead of doing 2 DFT's.


How to Show A OFDM Symbol's PSD in Matlab I'm trying to suppress the sidelobe of OFDM but having troubles showing its spectrum in Matlab. After IFFT, I get a time-domain signal, but in which way I can get its spectrum(or PSD)? I want to analyze the out band power emission, with rectangle window in the time domain the out of band power curve decrease slowly, so it may interference the adjacent channel user, to avoid that situation we have to apply some algorithms to suppress the out of band power emission caused by the sidelobe of OFDM signal. After apply the algorithm I have to plot some figures to demo the effect of my algorithm, so I need the whole PSD of OFDM signal and I also want the out of band part become more steep. I tried to use FFT and pwelch function of Matlab, but never get what I want. I doubt if I use them in right way. So I'm hoping someone here can help me. I know it's basic, but it really obsess for a long time. I analyze the baseband signal, so it should be after the IFFT operation.The figure that I want is below.

My matlab code:

  N = 64;
  over_sample_factor = 2;
  M = N*over_sample_factor;
  Mod = 16;
  symbol = 1;
  bitlength = N*log2(Mod)*symbol;
  itr_num = 1000;
  fft_len = 2*M;
  signal_freq = zeros(itr_num,fft_len);
  for itr = 1:itr_num
      bit_data = randi([0,1],bitlength,1);
      h = modem.qammod('M', Mod, 'SymbolOrder', 'Binary', 'InputType', 'Bit');
      zp_before = modulate(h,bit_data);
      after_zp = zeros(1,M);
      after_zp(1:N/2) = zp_before(N/2+1:N);
      after_zp(M-N/2+1:M) = zp_before(1:N/2);
      ofdm_symbol = ifft(after_zp);
      signal_freq(itr,:) = abs(fft(ofdm_symbol,fft_len)).^2;
  end
  PSD_mean = mean(signal_freq,1);
  plot(fftshift(10*log10(PSD_mean)));

The figure is like this

If you are looking at the sidelobes in particular then you will probably want to oversample, as Deve mentioned.  x2 oversampling should be sufficient, but you can do more if you want.  I generally plot my FFT's as follows to get the energy in a dB scale (Matlab syntax)-

plot(20*log10(abs(fft(dataVector))))

or, to center the FFT visually-

plot(fftshift(20*log10(abs(fft(dataVector)))))


confusion related to finding inverse Z-Transform using Complex Integral Method I am facing problems related to evaluation of inverse Z-Transform using Complex Integral Method;
Consider $X(z)=\frac{z}{z-2}  $ and $ROC: |z|>2$
then, $$x(n)= \frac{1}{2\pi j}\oint_c X(z) z^{n-1} \, dz$$
$$\implies x(n)= \frac{1}{2\pi j}\oint_c \frac{z^n}{z-2} \, dz \quad \dots (1)$$
$$\implies x(-1)=\frac{1}{2\pi j}\oint_c \frac{1}{z(z-2)} \, dz$$
$$=z\frac{1}{z(z-2)}|_{z=0} + (z-2)\frac{1}{z(z-2)}|_{z=2}$$
$$=-\frac{1}{2} + \frac{1}{2} =0$$
which is obvious as $x(n)=2^n u(n)$ (by inspection), $\implies x(-1)=0$
But for the same $X(z)$ , if $ROC: |z|<2$
then, from $eq(1)$ , we get:
$$x(1)=\frac{1}{2\pi j}\oint_c \frac{z}{z-2} \, dz \quad \dots (2)$$
$$=2 \quad \text{(from Cauchy's Integral Theorem)}$$
but $x(n)=-2^n u(-n-1)$ (by inspection), $\implies x(1)=0$
so,

How $c$ is defined in Complex Integral Method for a given $ROC$ ?
How can we evaluate $eq(2)$ ?

The closed contour $C$ must lie inside the region of convergence, so for the ROC $|z|<2$ you have no poles inside $C$ for $n\ge 0$, hence $x[n]=0$ for $n\ge 0$.


Why is peak detection in chromatography not completely automatic? When I started to write a software for chromatogram processing I naively thought I'd find some existing algorithms that do all the peak detection automatically and there won't be a need for any manual curation whatsoever. But it appeared that none of the major software in this field works this way. Users have to manually set up parameters (peak widths, thresholds) which work on some chromatograms but not others. So it's either an endless tuning of parameters or manual correction.
My question is philosophical: what keeps this problem unsolved? Do other fields (astronomy, particle physics) have similar problems interpreting measurements automatically? Can other approaches except DSP (e.g. Machine Learning) be useful here?

This is a problem we are also currently working on (added two recently published papers below), because we found it unresolved yet, and high-throughput experiments nowadays require more automated high-throughput analyses.
A first set of causes are DSP-like disturbances on the "interesting peak data": non-linearity (like saturation), added noise (with unknown distributions) and baselines (with unknown regularity). They can related to solute or matrix effects, type of chromatography (liquid, gas), columns, detectors, the complexity of the mixture, etc. We can include more or less visible processing performed inside the measurement tools: sampling (how any points per peaks?), smoothing, normalizations.
In gas chromatography there are models like Gaussian shapes, that can be explained by physical models (plates). However, separated peaks can exhibit distorsions such as fronting and tailing, or other asymmetries, that can vary throughout the chromatogram in uneven manners.
Even with "nicely"-sampled discrete Gaussians, very close peaks may be difficult to separate or distinguish (shouldering effects). More globally, deconvolution of sums of Gaussians remains difficult, because the direct problem is ill-posed: there are many ways to fit a single Gaussian with one of several Gaussians, and algorithms require additional optimization terms (constraints, regularization), that come with tuning parameters when combined with the data fidelity term (often a least-squares, but other choices are interesting).
As far as I know, such issues exists in other analytical chemistry modalities: NMR, infrared, mass spectrometry, Raman etc., each with their own issues. Variable peak shifting related to molecule interactions in NMR is a challenge for subsequent (blind) source separation algorithms.
There are works on using more machine learning techniques, yet it is likely that using physical or DSP models would help them.
A possible starting point on the difficulties of fitting: Dependence of Chromatogram Peak Areas Obtained by Curve-Fitting on the Choice of Peak Shape Function, 1997
We recently considered a couple of methods towards reconstruction or recovery of physico-chemical signals. They definitely are not the end of the game:
SPOQ ℓp-Over-ℓq Regularization for Sparse Signal Recovery applied to Mass Spectrometry, IEEE Transactions on Signal Processing, 2020 (Matlab code SPOQ toolbox)

Underdetermined or ill-posed inverse problems require additional
information for \ldd{d} sound solutions with tractable optimization
algorithms. Sparsity yields consequent heuristics to that matter, with
numerous applications in signal restoration, image recovery, or
machine learning. Since the ℓ0 count measure is barely tractable, many
statistical or learning approaches have invested in computable
proxies, such as the ℓ1 norm. However, the latter does not exhibit the
desirable property of scale invariance for sparse data. Extending the
SOOT Euclidean/Taxicab ℓ1-over-ℓ2 norm-ratio initially introduced for
blind deconvolution, we propose SPOQ, a family of smoothed
(approximately) scale-invariant penalty functions. It consists of a
Lipschitz-differentiable surrogate for ℓp-over-ℓq quasi-norm/norm
ratios with p∈]0,2[ and q≥2. This surrogate is embedded into a novel
majorize-minimize trust-region approach, generalizing the variable
metric forward-backward algorithm. For naturally sparse
mass-spectrometry signals, we show that SPOQ significantly outperforms
ℓ0, ℓ1, Cauchy, Welsch, SCAD and Celo penalties on several performance
measures. Guidelines on SPOQ hyperparameters tuning are also provided,
suggesting simple data-driven choices.

Sparse Signal Reconstruction for Nonlinear Models via Piecewise Rational Optimization, Signal Processing, 2021

We propose a method to reconstruct sparse signals degraded by a
nonlinear distortion and acquired at a limited sampling rate. Our
method formulates the reconstruction problem as a nonconvex
minimization of the sum of a data fitting term and a penalization
term. In contrast with most previous works which settle for
approximated local solutions, we seek for a global solution to the
obtained challenging nonconvex problem. Our global approach relies on
the so-called Lasserre relaxation of polynomial optimization. We here
specifically include in our approach the case of piecewise rational
functions, which makes it possible to address a wide class of
nonconvex exact and continuous relaxations of the ℓ0 penalization
function. Additionally, we study the complexity of the optimization
problem. It is shown how to use the structure of the problem to
lighten the computational burden efficiently. Finally, numerical
simulations illustrate the benefits of our method in terms of both
global optimality and signal reconstruction.


How to Create a Convolution Matrix with a Variable Condition Number (CN) I want to know the performance of a deconvolution algorithm with different CN, so I'm convolving my signal with different convolution matrices(different CNs) and then applying the deconvolution algorithm, then the error between the original and reconstructed signals is measured.

Is there a proper way to create the convolution matrix with variable CN.

Let's try to think of it intuitively.
Given the LPF what would make the inverse "Hard"?

The sections we won't be able to recover are where the LPF is 0, since there is no inverse for that we can multiply the result by.
Real world LPF won't reach "Real" zero as usually.
But what would make it hard to recover is where there are big ratios between the biggest and the smallest magnitude.
As close and fast your LPF goes to zero (And has magnitudes which are big) the harder the recovery will be.

Now just build analog LPF which this requirements at different scales, digitize it, create the Convolution Matrix and there you have it...

A practical example is given in my answer to 1D Deconvolution with Gaussian Kernel (MATLAB).


How to denoise a conversation between two people? I would like to know how to determine (using any programming language) whether a given audio file contains noise, such as of people talking in background, how to remove the noise, and how to subsequently compare the cleaned file with the original file.

The input to the program will be a sound file containing a conversation between two people and the program should be able to remove the noise from that file and should output only the conversation between the two people.

I think what you're attempting to do is a case of blind source separation (check Wikipedia). There are many techniques to solve this kind of problem, maybe you could check those on the Wikipedia and use the one that fits your needs best.

Particularly, last year one of my teachers showed us a demo of the ICA algorithm, and we got to play round with it in MATLAB. It seemed really effective, and there's a research group that even created this MATLAB toolbox for the ICA algorithm. Here you can check a simple demo (called cocktail party problem) of how ICA works and how effective it is.

Since you say you're new to DSP, I'm not sure if ICA or the other techniques would be too advanced for you, but maybe this will put you in the right direction. Good luck.


How to implement a $j$-level $M$-band wavelet transform of an image? I want to implement an $M$-band (or multi-band) wavelet transform, to be used for feature extraction on images.

Is there a built-in Matlab function available?
Do I have to implement is using filter coefficients?
Where can I get filter coefficients for $M$-band wavelets?

[EDIT: some code made available]  The  notion of $M$-band wavelet transform, $M\ge2$, generalizes the standard $2$-band wavelet. The theory is provides, for instance, in Theory Of Regular M-Band Wavelet Bases, P. Steffen et al., 1993. Remember that with $2$-band wavelets, one cannot obtain wavelets that are real-valued, symmetric, orthogonal and with finite support, except for the Haar wavelet (not very regular). For $M>2$, such designs can be obtained with sufficient regularity, and $M$-band wavelets are sometimes believed to be more efficient for image analysis.

Alas, regular $M$-band filterbank are generally harder to design, and $M$-band wavelets have been used more rarely. Resultingly, they are less implemented in standard software.

On $M$-band implementations: I am not aware of build-in Matlab functions. However, there exist $1$D versions, that you can built upon for your applications:

LTFAT (Large Time-Frequency Analysis Toolbox) : WFILT_CMBAND - Generates $$M$-Band cosine modulated wavelet filters, from On cosine-modulated wavelet orthonormal bases, 1995
M-band dual-tree wavelet: Shannon & Meyer wavelets, and Modulated Lapped Transforms for specific $M$s, and a $4$-band filter from Design of Efficient $M$-Band Coders with Linear-Phase and Perfect-Reconstruction Properties, 1995, Alkin and Caglar. They  come with their Hiilbert transform, and are detailed in Image Analysis Using a Dual-Tree M-Band Wavelet Transform, 2006; Noise Covariance Properties in Dual-Tree Wavelet Decompositions, 2007; A Nonlinear Stein Based Estimator for Multichannel Image Denoising, 2008. Below, 2D oriented $4$-band wavelets from Alkin and Caglar:

A $2$-D version is now embedded in the Matlab toolbox:  $M$-band $2$D dual-tree (Hilbert) wavelet multicomponent image denoising. Other sources to grab coefficients from:

-An algebraic construction of orthonormal $M$-band wavelets with perfect reconstruction

And one day, I shall rewrite my own codes to make them much cleaner.


Should we divide the power at 0Hz by two for Single-Sided Power Spectrum with FFT? To make my question as clear as possible I will go through an example. I want to plot one sided power spectrum of a signal. As an example my signal is sum of "50Hz sine with and amplitude of 1V" and "DC offset voltage of 1V". I create this signal and plot FFT of it where I define the power as amplitude square divided by two.
Here is my code:
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
import scipy.fftpack

f = 50 #signal freq
D = 1 #duration
fs = 800 #sampling freq
T = 1/f #signal period
N = int((D/(1/fs))+1) #number of smaples

t = np.linspace(0, D, num=N, endpoint=True) #time vector
dc = 1
y = np.sin(2*np.pi*f*t) + dc

plt.plot(t, y,'-b')
plt.plot(t, y,'.r')
plt.title('~ Sinusoid ~')
plt.xlabel('time [s]')
plt.ylabel('Voltage [V]')
plt.grid()
plt.show()

#FFT
plt.figure()
y = y
T = t[2]-t[1]
sampling_rate = 1/T
N = len(y)
yf = scipy.fftpack.fft(y)
xf = np.linspace(0.0, 1.0/(2.0*T), N//2)
amplitude = 2.0/N * np.abs(yf[:N//2])
pow = amplitude*amplitude/2

plt.plot(xf, pow,'b')
plt.grid()
plt.xlabel('Frequency [Hz]')
plt.ylabel('Power [W]')

Now the rms power for the sine component can be calculated as (amplitude square)/2 which is 0.5W and this is what we see at the plot above.
And for the 1V DC component, I would say the rms power is (amplitude square) which is 1W. But the plot shows twice of it namely 2W.
My question is: Should I divide the power by two at 0Hz at FFT plot in my code or am I interpreting something wrong?

The factor of "2" comes from the contribution of the negative frequencies, which doesn't apply to DC.
A better way to plot this is to NOT multiply with 2 but plot the entire FFT range from -400Hz to +400Hz. There you will see three components: 1 W at DC and 0.25W each at -50Hz and +50Hz which is exactly what's happening here.
If you want to plot power spectrum one sided, you need to multiply by two EXCEPT for DC and Nyquist, which have no negative frequency equivalent.


Correlation between amplitude and power In this link on page 218, just before the first paragraph ends, it says:
" We can interpret the power axes as one-half the square [...] of the modulus of the complex envelope."

In the same paragraph, the average input power is stated to be $P_\text{in} = A^2/2$. Can somebody please explain to me, why there is no resistance ?
The amplitude probably refers to a voltage signal. The RMS of a sinusoidal signal is $A/\sqrt(2)$. Since this is the square root of the mean, it makes sense to me, that the average input poewr is $A^2/2$. But still, I am confused about the lacking resistance. Especially, since in the second paragraph, it is stated, that the power axis will be graduated in units of $10\log(P_\text{in})$.

You're dealing with digital signals. Power hence doesn't have a physical unit, just as the amplitude of the signal doesn't have one – where it originated, i.e. the ADC, it might have had one, but usually, that unit just doesn't matter in DSP and the fact that the digital numbers are proportional with some constant factor to a physical signal is totally sufficient.

If you later need to actually map a digital power to a physical power, you'd just calibrate (ie. measure with a known source of amplitude or power) and calculate that constant factor.

TL;DR: This is DSP. We don't care about physical units, so (mentioning the) resistance is futile.


Academic problem solving in Computer Vision I am currently pursuing a Master degree of Engineering and I will have to write a thesis about solving a specific Machine Vision application. I have done some CV projects in the past but now I am tinkering about the methods of problem solving in the broad field of object detection. What I noticed is that it looks like there are about two kind of people trying to solve CV problems and asking about them on a SE site.

The first kind of people needs to solve a specific problem (like, finding faces or horses or else in a picture) and is asking questions like "Oh, I need to find an algorithm to detect face" or "How can I detect cars with SIFT and OPENCV".

The second kind of people is writing academic papers on some new algorithms (like SIFT, ORB, etc.). They might have questions like "How to best compute the Hamming distance" or else.

Both kind of people and questions are relevant and fine to me (if my opinion matters). I wonder though if there is somewhere I could ask questions or share thoughts about more structured methods of computer vision problem solving. I mean, I would like to know a place where I could ask questions such as:

"What are the general steps of solving an object detection (or recognition) in a 'graduate' manner, and being able to demonstrate that all the steps used where necessary and to a high confidence the best method was chosen?"

I wanted to post on Signal Processing Meta but you need 5 rep in order to...

BTW, the question asked is a real one that I would like to ask here, simply it might not be the place or respect the guidelines of this SE site. In the same way, I feel this might be too CV related to ask on Academia SE...

Edit
What I mean by "a graduate manner" is the rigorous approach to problem solving that one should use in the course of graduate studies. What I mean is that sometime (especially when you work in a company vs working in academia), you must rapidly get a solution that work to solve a problem, instead of doing an entire review of the possible methods and doing benchmarks and etc... (This is a generalization and I don't want to imply that working in a company makes you use somewhat botched methodology, see it more like a "Perfect is the enemy of good".)

With all due respect, your question is rather vague...  For example, it is not clear to me what you mean by "structured methods of computer vision problem solving", or what you mean by a "graduate manner" of object detection or recognition.

If you are doing a Master's thesis in computer vision then presumably you have an adviser, who is supposed to be knowledgeable in this field, and who is supposed to help you define your problem, and point you to the relevant papers.

By the way, defining and constraining your problem is the key here. "general steps of solving an object detection (or recognition)" is not a very well defined problem.  Are you trying to detect a specific object, such as a particular book or your favorite coffee mug? Or are you trying to detect members of broad categories of objects, such as faces, cars, people, or chairs?

Does your method need to be invariant to in-plane rotation?  The answer is probably "no" for faces, which are rarely upside down, but "yes" for detecting cars in aerial images.  There are many more questions like this.

So, my point is that object detection, let alone recognition, is a broad area, and it is best to talk to your adviser about scoping it down.


SIS Particle Filter Python Implementation Lets say I have some time series data which I generated like this:

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from numpy.random import uniform

# Length of time series
N=400

# Gaussian random numbers as an excitation signal
ex = np.random.randn(N)

# Second order AR Process with coefficients slowly changing in time
a0 = np.array([1.2, -0.4])
A = np.zeros((N,2))
omega, alpha = N/2, 0.1

for n in range(N):
    A[n,0] = a0[0] + alpha * np.cos(2*np.pi*n/N)
    A[n,1] = a0[1] + alpha * np.sin(np.pi*n/N)

S = ex.copy()
for n in range(2,N):
    x = np.array([S[n-1], S[n-2]])
    S[n] = np.dot(x, A[n,:]) + ex[n]

fig,ax = plt.subplots(nrows=2, ncols=2, figsize=(9,4))
plt.tight_layout()
ax[1,0].plot(range(N), A[:,0])
ax[1,0].grid(True)
ax[1,0].set_title("Coefficient a0", color='m')
ax[1,1].plot(range(N), A[:,1], color='m')
ax[1,1].grid(True)
ax[1,1].set_title("Coefficient a1", color='m')
ax[0,0].plot(range(N), ex)
ax[0,0].grid(True)
ax[0,0].set_title("Random Excitation Signal")
ax[0,1].plot(range(N), S, color='m')
ax[0,1].grid(True)
ax[0,1].set_title("Time Varying Autoregressive Process")

How do I implement a Sequential Importance Sampling Particle Filter on this data? I understand all the theory and the maths behind it and even saw some examples of implementation on different problems but none of them work with similar data. For example, https://github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python/blob/master/12-Particle-Filters.ipynb deals with robot movements and landmarks and I just want to apply the algorithm to a simple time series data...Can anyone help me with the implementation? I'm especially confused as to how to do the sampling part...

Particle filter is a Bayesian filter. Any Bayesian filter requires process and measurement models so you also need to define them.
If you will use bootstrap particle filter, you just create initial samples with initial parameters from Gaussian distribution. Then, you propagate the particles using process model. This is your so called "prediction" step.
After obtaining your measurements for each time step, you calculate the weights. Using those weights you perform resampling.


Weird results out of fourier transform of unknown signal, what am I doing wrong? It has been a long time since I studied engineering, please forgive my mistakes.

I got data sampled from a digital accelerometer. It is sampled at 20 hz, then collapsed into 1 second epochs. The data is then filtered to remove noise.  The system outputs 0 at rest, a lot more during activity.

I'm looking at the fourier transform of the signal, 1024 samples, sampled once per second. I'm  seeing this:

Here's a real-imaginary plot of the real-imaginary results of the Fourier transform.

The signal appears to be a helix in frequency domain.The system appears to oscillate around some limit cycle within frequency domain.

Is this possible?
Is Fourier transform completely wrong for the digital signal?
 Should using only the DFT for this kind of analysis?
Is my window size wrong ?
 Am I seeing aliasing due to incorrect window size?

It has been a few years since my last engineering class, and I'd really appreciate any help on the subject. I got 3 books on engineering, signal analysis and transforms, but it takes time for this kind of knowledge to come back to me.

My guess: nothing to see here. Move along.

You mention that "the signal appears to be a helix in [the] frequency domain". I assert that what this helix really is, is a complex exponential. And there is nothing wrong with that, because complex exponentials are what you get when you time-delay a signal: $f(x-a)\Leftrightarrow \hat{f}(\xi)e^{-2\pi i a \xi}$. The helix pattern simply reflects that most of the power in the signal you are acquiring isn't located at $T=0$. (And it shouldn't be there anyway!)

Eyeballing, the period of oscillation in your frequency plot appears to be ~28 samples. At a sampling rate of 20hz, if what I'm saying is true, then most of the energy in your time-domain signal should be centered roughly around $t\approx\frac{28}{20}=1.4$ seconds. Was I close?

The real problem here is likely to be your plots. It seems like you're plotting raw real/imag DFT outputs. Don't do that. Instead, compute magnitude/phase from real/imag values, and plot that.


Can I use SINAD for THD+N measurements? Are they inverse of each other? I'm looking for a clarification, which I haven't really been able to find. It's probably obvious but I'd like to be absolutely sure of what I'm calculating.

I noticed that the formula that calculates THD is:

$THD = \frac{\sqrt{H^2_2 + H^2_3 + H^2_4 + \dots + H^2_n}}{H_1}$

This is seen as the power ratio of the harmonics to the fundamental.

The sinad function in MATLAB calculates the power ratio of the fundamental to the harmonics and noise floor.

$SINAD = \frac{P_{fundamental}}{P_{noise+distortion}}$

Wouldn't THD+N be the inverse of SINAD, since for calculating THD+N we essentially need:

$\frac{P_{noise+distortion}}{P_{fundamental}}$

Both the measurements seem to ignore DC. I use this MATLAB tutorial page as a reference: https://www.mathworks.com/help/signal/ug/distortion-measurements.html

I guess I'm wondering this because I often see SINAD and THD+N listed as separate measurements, but I don't see how they are different.

So I guess I should have looked around a little better. I've found that SINAD and THD+N are reciprocals.

These references were useful:

https://www.ap.com/technical-library/more-about-thdn-and-thd/

http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.476.1903&rep=rep1&type=pdf


Identify which of the three signals is closest to a sinusoidal curve I have some acceleration values from a sensor. These values are reported as x, y and z components of the acceleration. As can be seen in the attached figure, it is clear that the z signal has a clear sinusoidal curve present. Is there any algorithm to quantify the "sinusoidalness" of a curve. I am trying to write an algorithm that would identify the z signal of the attached image to be the most sinusoidal.

Another example is shown in the figure below: For this image the x-axis contains the < most sinusoidal > signal and the algorithm must be able to identify the x-axis as the correct axis.

Edit: The variation could indeed be greater in any of the axes without the most sinusoidal signal. The original problem being solved here is to count the number of breaths in the accelerometer signal. To do this first the axis containing the breathing signal must be identified. From the image it is obvious that for the first image the z-axis has the breathing data and for the second image the x-axis has the breathing data

Is it sufficient to identify the «most sinuoidal» of those 3, or would you also want linear projections of those (consistent with a IMU sensor tilted vs the plane of motion)?
A simple solution might be to do a windowed fft and pick the direction where the «crest factor» of the fft magnitude was largest (best explained by a single sinoid).
Edit:
It appears thah your «sine» is somewhat variable. If the above approach fails, then perhaps removing the DC component («highpass filter») followed by removing the high frequency edges and impulses (lowpass filter) for an effective tuned bandpass filter, followed by squaring and summing the samples could be used.


Digital Image Processing I am trying to extract only urban aeas from satellite imagery Using Matlab. I have applied some edge detectors like canny and also used morphological functions and I am able to extract all the boundaries ,but not able to seperate them from various urban aeas like roads, trees and waterbodies.

Please help i can share the code if any one wants

I think this code would be helpful, as displaying intensity part gives you the spatial characteristics of an image.

%Write a MATLAB function which converts an RGB to HSI. Display the Hue
%image, Saturation image and the Intensity image. Histogram equalize the H,
%S and I images and the HSI image. Then convert the HSI image back to RGB.

function rgbtohsi(x)
F=imread('C:\Documents and Settings\admin\Desktop\TM\tanya\img3.tif ');
F=im2double(F);
r=F(:,:,1);
g=F(:,:,2);
b=F(:,:,3);
th=acos((0.5*((r-g)+(r-b)))./((sqrt((r-g).^2+(r-b).*(g-b)))+eps));
H=th;
H(b>g)=2*pi-H(b>g);
H=H/(2*pi);
S=1-3.*(min(min(r,g),b))./(r+g+b+eps);
I=(r+g+b)/3;
hsi=cat(3,H,S,I);
HE=H*2*pi;
HE=histeq(HE);
HE=HE/(2*pi);
SE=histeq(S);
IE=histeq(I);
choice=input('1:RGB to HSI\n2:Display Hue, Saturation and Intensity Images\n3:HSI to
    RGB\n4:Hue-Equalization\n5:Saturation-Equalization\n6:Intensity-
    Equalization\n7:HSI-   Equalization\n Enter your choice :');
switch choice
case 1
    figure,imshow(F),title('RGB Image');
    figure, imshow(hsi),title('HSI Image');
case 2
    figure,imshow(F),title('RGB Image');
    figure, imshow(H),title('Hue Image');
    figure, imshow(S),title('Saturation Image');
    figure, imshow(I),title('Intensity Image');
case 3
    C=hsitorgb(hsi);
    figure,imshow(hsi),title('HSI Image');
    figure, imshow(C),title('RGB Image');
case 4
    RV=cat(3,HE,S,I);
    C=hsitorgb(RV);
    figure,imshow(hsi),title('HSI Image');
    figure,imshow(F),title('RGB Image');
    figure, imshow(C),title('RGB Image-Hue Equalized');
case 5
    RV=cat(3,H,SE,I);
    C=hsitorgb(RV);
    figure,imshow(hsi),title('HSI Image');
    figure,imshow(F),title('RGB Image');
    figure, imshow(C),title('RGB Image-Saturation Equalized');
case 6
    RV=cat(3,H,S,IE);
    C=hsitorgb(RV);
    figure,imshow(hsi),title('HSI Image');
    figure,imshow(F),title('RGB Image');
    figure, imshow(C),title('RGB Image-Intensity Equalized');
case 7
    RV=cat(3,HE,SE,IE);
    C=hsitorgb(RV);
    figure,imshow(hsi),title('HSI Image');
    figure,imshow(F),title('RGB Image');
    figure, imshow(C),title('RGB Image-HSI Equalized');
    F1=C;
    otherwise
    display('Wrong choice');
end
end
    function C=hsitorgb(hsi)
    HV=hsi(:,:,1)*2*pi;
    SV=hsi(:,:,2);
    IV=hsi(:,:,3);
    R=zeros(size(HV));
    G=zeros(size(HV));
    B=zeros(size(HV));
%RG Sector
    id=find((0<=HV)& (HV<2*pi/3));
    B(id)=IV(id).*(1-SV(id));
    R(id)=IV(id).*(1+SV(id).*cos(HV(id))./cos(pi/3-HV(id)));
    G(id)=3*IV(id)-(R(id)+B(id));

%BG Sector
    id=find((2*pi/3<=HV)& (HV<4*pi/3));
    R(id)=IV(id).*(1-SV(id));
    G(id)=IV(id).*(1+SV(id).*cos(HV(id)-2*pi/3)./cos(pi-HV(id)));
    B(id)=3*IV(id)-(R(id)+G(id));
%BR Sector
    id=find((4*pi/3<=HV)& (HV<2*pi));
    G(id)=IV(id).*(1-SV(id));
    B(id)=IV(id).*(1+SV(id).*cos(HV(id)-4*pi/3)./cos(5*pi/3-HV(id)));
    R(id)=3*IV(id)-(G(id)+B(id));
    C=cat(3,R,G,B);
    C=max(min(C,1),0);
end