Correlation between a signal segment and its stretched replica In my application the measurements are affected by temperature and the signal is stretched over time, though preserving relatively similar structure.

I want to find incremental stretching between pairs of consecutive signals. Having a pair of signals, I divide the first into N segments (windows) and for each window I want to find it's shift and new length in the next signal.

I tried to estimate these incremental shifts with peak detection, but got low accuracy estimates, probably due to the stretching effects on the waveform structure.

In addition, I was advised to try Dynamic Time Warping algorithm, but it manipulates both signals to find the minimal Euclidean distance. In my case I would want only the first segment to be warped until it correlates perfectly with some waveform in the 2nd signal. Couldn't find the way to do so with DTW.

Say you got two similar signals when one of them is stretched and you want to know the straching ratio...
Seems like you can
1. guess the straching ratio
2. resample the short signal with a resampling ratio equal to your guessed straching ratio.
3. Find the max of the convolution between the signals.
4. Optimize the estimation of the straching ratio using a numerical method.


detect to rising, stable and falling point in non-smooth rectangular wave I am working on basic signal processing problems in MATLAB. I have found a signal from the internet (i don't remember the site exactly). The data is organized in column wise. 1st and 2nd column is data and 3rd column is a rectangular wave. The rectangular wave is not very smooth. So I decided to detect four points as shown in the figure. I spend 2 days to find the exact solution but i couldn't find any solution to the detect those four points.
I tried few things:
1) finding the difference between the points to get 1,2,3 and 4 points.
but it doesn't work.
2) Median Filter to smooth the wave.

Anyone, please give some suggestion to detect these points in MATLAB

The usual approach to change detection is the CUSUM algorithm.

I've done an implementation that just addresses the level (mean) change issue. It's included (in R) below.

The black line is the noise-free data, the red line is the noisy data and the blue bars are the detected breaks (for this realization).

This just addresses the level change; to address the three zones: low level, changing level, and high level, you'll need to figure out the way to estimate the mean in the changing level (perhaps assume a rise time and fix the two levels?).

The code below is based on this document.

A slight change might get you the first and second changes easily: instead of working with the data, work with the difference of the data.  If I do that, then I get:

but I've had to add:

sigma <- 0.001
data_before_diff <- noiseless_data + rnorm(length(noiseless_data),0,sigma)
data <- diff(data_before_diff)
mu2diff <- 0.04
thresh <- 1

and then later:

 s[k] = (mu2diff - mu1)/sigma*(data[k] - (mu2diff + mu1)/2)

and

    tmp <- mu1
    mu1 <- mu2diff
    mu2diff <- tmp

To find the changes at the end of the signal, you'll probably need to change mu2diff to -mu2diff and check for that versus 0.

Note I've had to greatly reduce the noise variance and thus the threshold to get this to work.

R Code Below

# 30039
# CUSUM

N <- 30

mu1 <- 0
mu2 <- 1
noiseless_data <- c(rep(mu1,N), seq(mu1,mu2,1/N), rep(mu2,N*4), seq(mu2,mu1,-1/N), rep(mu1,N))

sigma <- 0.1
data <- noiseless_data + rnorm(length(noiseless_data),0,sigma)

thresh <- 100

breaks <- rep(0,N)
num_breaks <- 0

s <- rep(0,length(data))
capS <- rep(0,length(data))
G <- rep(0,length(data))
for (k in 1:length(data))
{
  s[k] = (mu2 - mu1)/sigma*(data[k] - (mu2+mu1)/2)
  if (k==1)
  {
    capS[k] = s[k]
  }
  else
  {
    capS[k] = capS[k-1] + s[k]
  }

  G[k] <- max(0,capS[k] - min(capS[1:k]))

  if (abs(G[k]) > thresh)
  {
    #capS[k] <- 0
    num_breaks <- num_breaks + 1
    breaks[num_breaks] <- which.min(capS[1:k])
    tmp <- mu1
    mu1 <- mu2
    mu2 <- tmp
  }
}

first_break <- min(breaks[1:num_breaks])
last_break <- max(breaks[1:num_breaks])

plot(data,col="red", type="l")
lines(noiseless_data)
lines(c(first_break,first_break),c(0,1),col="blue", lwd=10)
lines(c(last_break,last_break),c(0,1),col="blue",lwd=10)


How detect and filter noise from signal This is the spectrum of the signal:

And this is the signal:

On the FFT spectrum you can notice the high amplitude peaks and I think these frequencies are the main ones in the signal:

I think that the "low" harmonic must be clean, but I have the problem that low amplitude components are placed in every frequency interval.  The main problem is that I haven't got any model of the signal for a hypothesis about noise, and as result, I must detect noise components.  After this step, I must design a filter or use another method of filtering for denoising the signal.  Which methods can I apply for (1) detecting noise and (2) filtering?

Just "eyeballing" your signal, it looks like the interesting things are all down below 1Hz.  The peaks in your signal browser diagramm seem to be spaced at every 5 to 10 Seconds, which would be 0.1 to 0.2 Hz.  Your spectrum plot goes from zero to 25 Hz, and doesn't show anything really interesting.

Try zooming in on the spectrum and see what lies below 1 Hz.  Yo may well find that the signal itself lives down in the very low frequencies, and that separating the noise from the signal is as simple as using a low pass filter.


What is an intuitive explanation of cepstrum for a 12-year-old? In Feynman's spirit of "If you truly understand a concept, you should be able to explain it clearly" and the explanation should make solid sense.

There are two main challenges:

The axis of "Quefrency" is not a straightforward/natural concept, how to describe and make sense of what it is?
Why is the cepstral view an efficieny way of "elevating" the main information contained in the signal? We can see that the vector is more succinct, but where does this efficiency come from?

“efficiency” comes from a-priori assumptions.  One assumption that makes a Cepstrum useful is that the signal contains something roughly periodic including a set of harmonics or overtones, possibly  inside some sort of frequency spectrum envelope.

So a Cepstrum is mostly just an FFT of an FFT.  with a few mods.  On mod is to take the log magnitude after the first FFT to keep big overtones from drowning out the little ones.

So a DFT (slow FFT) correlates periodic happenings against some sinusoidal basis vectors.  Unless the periodic time domain  happenings are a single pure sinusoid, those periodic happenings will often correlate well against multiple equally spaced overtone frequency DFT basis vectors.

Thus you end up where you started, with a bunch of roughly equally spaced happenings, but now in the frequency domain instead of the time domain.  The second FFT will again correlate these against sinusoidal basis vectors related to that spectral periodicity, which for audio signals turns out to be an interesting proxy for (overtone rich) pitch.  The spectral envelope, since it was forced to be positive (by the log magnitude), ends up surrounding DC, instead of becoming sidebands of the spectral peaks as with time domain modulation.  This helps separate envelope information (formant, etc.) from pitch information in the Cepstrum result.

A Cepstrum also works for a single pure sinusoid, as the two spectral lines (positive and negative frequency peaks) correlate against a negative cosine of the same width as the line separation.

Note a required trick: One often has to make zero and tiny FFT spectrum magnitudes bigger before the log  to avoid the log from producing gigantic negative spikes or NaNs.


Determinant of Hessian approximation (SURF) I have a question regarding formula in SURF article by Bay et al.

Theory

Given a point $p=(x,y)$ in an image $I$, the Hessian matrix $\mathcal{H}$ in $x$ at scale $\sigma$ is defined as follows

$$
\mathcal{H}(p, \sigma) = \begin{bmatrix}L_{xx}(p, \sigma) & L_{xy}(p, \sigma)\\L_{xy}(p, \sigma)& L_{yy}(p, \sigma)\end{bmatrix},
$$

where $L_{xx}(p,\sigma)$ is the convolution of the Gaussian second
order derivative $\frac{\partial^2}{\partial x^2}g(\sigma)$ with the image $I$ in point p, and similarly for $L_{xy}(p,\sigma), L_{yy}(p,\sigma)$.

The $9\times9$ box filters in Fig. 2 are approximations of a
Gaussian with $\sigma=1.2$ and represent the lowest scale (i.e.
highest spatial resolution) for computing the blob response
maps. We will denote them by $D_{xx}$, $D_{yy}$, and $D_{xy}$. The
weights applied to the rectangular regions are kept simple
for computational efficiency. This yields
$$
\det(\mathcal{H}_{approx}) = D_{xx}D_{yy} - (wD_{xy})^2.
$$

The relative weight $w$ of the filter responses is used to balance
the expression for the Hessian’s determinant. This is
needed for the energy conservation between the Gaussian
kernels and the approximated Gaussian kernels,
$$
w = \frac{\big\lvert L_{xy}(1.2)\big\rvert_F\big\lvert D_{yy}(9)\big\rvert_F}{\big\lvert L_{yy}(1.2)\big\rvert_F\big\lvert D_{xy}(9)\big\rvert_F} = 0.912\ldots \approx 0.9,
$$

where $|X|_F$ is the Frobenius norm. Notice that for theoretical
correctness, the weighting changes depending on the
scale. In practice, we keep this factor constant, as this did
not have a significant impact on the results in our
experiments.

Questions

Main question is how this formula for $w$ is obtained?
And my another concern is to why $w$ is within the brackets with $D_{xy}$. If it's just a balancing coefficient wouldn't it make more sense to write like this $wD_{xy}^2$? So I guess there is some motivation behind it.

My general purpose is to generalize this method for 3D, so it would be great if anyone could share some thoughts about balancing coefficients in that case or some useful links containing relevant information.

I think $w$ is a factor that makes the ratio of  $D_{xy}$ to $D_{yy}$ the same as $L_{xy}$ to $L_{yy}$. This is so the value of the determinant for the simplified kernel roughly matches that of the continuous version.
You could write it that way. But if you leave it as $wD_{xy}$ then if used in other expressions you don't have to muck around with adding square roots. e.g. $(wD_{xx})^3$ expressed your suggested way would be $\sqrt{w}^3D_{xx}^3$


Can the FFT representation of a signal be modified w/o creating glitches? I've read that modifying the bin values in the frequency domain after a FFT poses problems that lead to glitches when one takes the signal back to time domain using iFFT. This is caused iirc because the modified data does not align "continuously" anymore to the bins (and discontinuities $\implies$ artefacts). Because of this it's e.g. not recommended to perform ("ideal filter") equalisation in the frequency domain.

But since there are e.g. interpolation methods,

Are there ways for doing processing such as shuffling or removing bins in the frequency domain such that one can retain sound quality? E.g. by "synthesising" new data using interpolation methods.

In the context of audio something seems to be feasible, since there exists plug-ins such as this that seem to do filtering in the frequency domain:

The way you do it is you use a window that tapers to 0 at the ends, both before and after the FFT (analysis and re-synthesis windows), along with 50% overlap (needs to be a COLA window), to fade out the glitches and keep the signal continuous.


Fractional Delay Filters and Cutoff Frequences I am trying to implement an interpolator for arbitrary sampling rate conversion of a one-dimensional signal (a fractional delay filter). I am aware of the fact that interpolation is in general a problem of low-pass filtering, and impulse response and passband widths are inversely proportionate (for the same kernel shape). However, I'm a bit uncertain of how precisely to apply this to fractional delay filters for different interpolation rates.

As I understand it, the frequency cutoff for the interpolation kernel should be at pi*low_rate/high_rate radians. But for what class of interpolation filters is this the case? E.g. for cubic B-spline interpolation (a non-Nyquist(1) piecewise polynomial, where Nyquist(1) means the filter's impulse response is delta when the source and target sampling rates are identical) for supersampling the kernel is based purely on the fractional position, and does not have any kind of provision for a cutoff, despite it being possible to compute an impulse response for the interpolator (which would seem to make the kernel frequency-dependent by the impulse response/passband duality mentioned previously). What about for other classes of interpolators such as linear (Nyquist(1), polynomial), sinc (Nyquist(1), non-polynomial) or Gaussian (non-Nyquist(1), non-polynomial)?

Additionally, non-Nyquist(1) interpolators can by converted to Nyquist(1) with the use of a fixed prefilter, e.g. as described in Variable Sample Rate Conversion Techniques for the Advanced Receiver. If in fact a given interpolator must have its impulse response dilated/contracted based on cutoff frequency, this means there are now two stages for these non-Nyquist(1) interpolators: the fixed prefilter and the potentially variable interpolation. At what stage(s) would the correction for cutoff frequency be applied in such a system, and would it follow the same formula stated above? I would have expected frequency correction to be necessary in the prefilter, but the paper doesn't seem to indicate the use of frequency-dependent prefilters with cubic B-spline interpolation despite that interpolator, as already stated, having no accommodation for cutoff frequency in the core algorithm.

Okay, I think I've learned a sufficient amount about interpolation theory to mostly answer this. Interpolation can be viewed in two equivalent ways. The most obvious is, for each point you are trying to compute the interpolated value for, the multiplication/sum of the source samples with the interpolation kernel centered on the interpolation point. Alternately, you could view the signal itself as the convolution of a set of "hidden" coefficients with the very same convolution kernel; in this case to interpolate between samples you multiply/sum each source coefficient by the point on that sample's corresponding interpolation kernel where the interpolation point falls. The former model is what you actually use for interpolating a signal, but the second makes obvious some things that are less clear in the first.

For upsampling, the equivalent low-pass filter cutoff will always be at the Nyquist rate of the source signal (at least for an ideal low-pass filter), meaning that the interpolation kernel will be fixed regardless of upsampling ratio. Furthermore, it means that the coefficients for a given interpolation point will always have fixed distance on the interpolation kernel curve, which is valuable information for optimization (e.g. if you quantize your frequency ratio to the same resolution as the kernel LUT, you can arrange the kernel LUT such that all coefficients for a single interpolated sample can be fetched at once, rather than requiring numerous array accesses per interpolated sample). And as a consequence of these, the prefilter (if needed) will be fixed as well.

For downsampling, the cutoff frequency must be the Nyquist rate for the target sampling rate, meaning the interpolation kernel as well as the prefilter must be dilated proportionately when applied to the source samples. I am not entirely certain, but I think this would mean interpolation methods whose kernels cannot be dilated (e.g. matrix-based B-spline interpolation) cannot directly be used for downsampling, in which case you would have to use such an algorithm to upsample to some integer multiple of the target frequency then low-pass filter and decimate (obviously simultaneous polyphase low-pass filter and decimate would be an optimization).


How to determine the parameter of a Complementary Filter? I know that the Complementary Filter has the functions of both LPF and HPF. But I think my understanding on the principal behind it is still unclear.

I am quite new on digital signal processing, and maybe some very fundamental explanations will help a lot.

Say I have a Complementary Filter as follows:

$$y =a\times y+(1-a)\times x$$

Then my parameter $a$ may be calculated by $$a=\frac{time\ constant}{time\ constant+sample\ period}$$
where the $sample\ period$ is simply the reciprocal of the $sampling\ frequency$.

The $time\ constant$ seems to be at my own choice.

My Questions:

What is the theory behind this calculation?
How do we choose the $time\ constant$ properly?

This is the difference equation for a low pass filter.  Five time constants (5 * 1 time constant) is the time it will take to for the output, to reach 99.33% of the value of the input, from when the input changes from 0 to its final value, and stays there (a step response).  This comes from 1-e^(-5).  As you are specifying time constants, it implies you want a time domain definition for use as a moving average filter, so specify the period of averaging you would like, assume that as 5 time constants and divide by 5 five to get one time constant for use in the equation you have quoted above.


Are there analogues to orthogonal transformations in non-orientable surfaces? I am working on extremely large, symmetric matrices of counts, and attempting to identify patterns/shapes within them.  Wavelets are a popular tool in image processing, and have some nice statistical properties when applies to discrete random variable.

The problem I face is this:  2D wavelet transforms on the diagonal contain redundant values, making the statistical interpretation of their coefficients difficult.

HOWEVER, i noticed that the lower triangular portion of my dataset (i.e. the dataset without any redundant values) has a natural analogue as a Möbius strip!

The steps taken to remove the redundant entries in my matrix are illustrated in the image below:

Remove the upper right half of the matrix.
Slice the image diagonally through the middle.
Twist and flip the lower right portion and "glue" it to the upper left portion so that the dashed black lines are aligned and pointing in the same direction.
Give the square a half twist and "glue" the yellow sides so that the arrows align in the same direction.

If you're incredulous that these steps are justified, see this video about redundant ordered pairs.

This is all well an good, but is there a natural extension of wavelets (or ANY orthogonal transformation) onto a non-orientable surface?  If so, how does one interpret a multi-resolution analysis?

My intuition is that horizontal detail wavelets will carry over, but the non-orientable property will complicate vertical and diagonal detail wavelets.

I am interested in sources that discuss the topic.  Web searches haven't turned up anything. I have emailed the topology faculty at my school, and have yet to hear back.

Bear with me, as this is really only tangential¹ to my area of expertise, but:

Maybe the question is not how to extend wavelets to work on non-orientable surfaces, but to interpret said surface as something where wavelets natively work well.

Enter differential geometry (differential topology? I don't know.):

As you've notice: if you only consider a part of your moebius strip, let's say, everything but the original "yellow line", then it looks perfectly much like structure that you can work with: a (discrete) two-dimensional vector space. I presume you have some explicitly vector-space notation of what your matrix represents!

In fact, the idea here would be to treat your surface as a vector bundle:
You take a circle (which is, if you think about it, a one-dimensional vector space itself, with some kind of "modulo" characteristics). To every point on that circle, you assign yet another vector space: a one dimensional line, tangential to that circle.

Now, take a part of your circle (an arc), and find the product space of that with your vector bundle: That's something that should look topologically like your original vector space.

This way of looking at your moebius strip should allow you to use your transforms like you used to use it on the original matrix.

¹ hrhr, a pun.


Alternating sign after Discrete Fourier transform Let me start by saying that this question is motivated by data analysis that I am doing. Its part of a much larger set of signal processing operations, but in order to keep the question contained I'll talk only about what (to my best knowledge) is relevant.

The scenario is as follows: we start from some discrete dataset $u$ of size $N$ in the time domain, $u ={u_1,u_2,\dots,u_N}$. From observation one finds that this dataset has some special properties: every odd value is purely real, and every even value is purely imaginary. If I subsequently compute the discrete Fourier transform of this dataset given by $v ={v_1,v_2,\dots,v_N}$ where
\begin{equation}
v_s=\frac{1}{\sqrt{N}}\sum_{r=1}^{N}u_r\exp{\left[2\pi i (r-1)(s-1)/n\right]}
\end{equation}

I find that I obtain a curious result: the outcome $v$ is almost purely real, and the even and odd values alternate in sign: $\mathrm{sgn}{(v_s)}=-\mathrm{sgn}{(v_{s+1})}$. With this first part I mean that the real data is $\mathcal{O}(10^{-4})$ and the imaginary part is $\mathcal{O}(10^{-12})$; I'd say is purely from numerical errors.

However, it is this second aspect that is not obvious to me; why do the real values alternate in sign? For those who appreciate some visual aspects, the real part of the Fourier transform ($\mathrm{Re}(v)$) looks like this
. Note that the first quarter is a mirror image of the second, and the third of the fourth; this is a different point that is not entirely relevant.

Note that an answer could also definitely be that, given the above, there is no explanation. I am not posting the entire dataset as its rather large and makes the question messy, but perhaps I am missing a crucial ingredient. In that case I will edit it to make it more complete.

Let's suppose that your original data, $u_n$, has been windowed (which is something you might typically do, in order to 'tidy up' the FFT). In this case, most of the energy in your input data sequence is clustered around the region near $n=N/2$. To be more specific, I expect that your input data is actually centred on $n=N/2+1$.

Now, look at your expression for the FFT, and imagine that all the $u_n$ are zero, except $u_{N/2+1}=1$. This means that:

$v_s=\frac{1}{N}\sum_{r=1}^N u_r \exp[2\pi i(r−1)(s−1)/n] = \frac{1}{N} \exp[2\pi i(N/2)(s-1)] = \frac{1}{N} (-1)^{s-1}$

In other words, the sign-flipping behaviour of your FFT output is generally related to the fact that the energy of you input sequence, $u_n$, is clustered around time $N/2$.

Ideally, you'd like to have your energy clustered around time zero, and in order to do this, you could shift your signal $N/2$ samples to the left (this should be a circular shift, so that the data wraps around in  the buffer). In practice, you don't bother with the circular shift, and instead you simply swap the sign of all the odd FFT bins, which does exactly the same thing as shifting the data in the time domain.

Also, if your FFT output is real, this implies that your time-domain input is Hermitian. This means that there is a particular symmetry in your input data, $u_n$, such that the real part of the input satisfies:

$Re(u_{N/2+1+k}) = Re(u_{N/2+1-k})$ for all $k$ $\quad\leftarrow$ This means that the real part is even

and the odd part satisfies:

$Im(u_{N/2+1+k}) = -Im(u_{N/2+1-k})$ for all $k$ $\quad\leftarrow$ This means that the imaginary part is odd

If you look at your data (the $u_n$ values) you should see this property. You will find that $u_{N/2+1-2}$ (which is real, according to your description) will be equal to $u_{N/2+1+2}$. Furthermore, $u_{N/2+1-1}$ (which is imaginary, according to your description) will be equal to $-u_{N/2+1+1}$.


FIR Filter - Transfer Function I am dealing with the transfer function of an FIR filter:

  $$H(z) = (1-0.5z^{-1})(1-2z^{-1})$$

I am having trouble determining which type of Linear-phase FIR filter it is, Type 1-4. I believe it may be Type 2, but not certain.

Also I need to determine the frequency response of the filter. Is it just as simple as replacing the $z^{-1}$ terms with $e^{-j\omega}$ in order to get the frequency response? If so this results in frequency response of

  $$H(e^{j\omega})=e^{-j\omega}(-2.5 + 2\cos(\omega))$$

And I am not sure how to take the magnitude or phase of that frequency response, any help would be greatly appreciated.

HINTS:

Just multiply the terms to get $H(z)$ in the form $$H(z)=h[0]+h[1]z^{-1}+h[2]z^{-2}$$ If you know the four types then you should immediately see which type it is.
Your frequency response looks good. The term in parentheses is real-valued. Does it change sign? You can directly write down magnitude and phase from that representation (just watch out for a negative sign, and what that means for the phase).


Two-GSMK channels demodulation: 3dB difference in BER I am trying to demodulate individually two GSMK channels separeted by 50kHz (one channel = 25kHz).
I am working on a limited ressources hardware so I decided to make some optimization on the signal processing chain.

Here is the architecture: after sampling, my two channels are aliased at 4.6MHz and at 4.65MHz.

The architecture for one channel is as follow: a DDS (Direct Digital Synthesis) and a mixer bring my signal in baseband. I then have a downconversion, filtering, the decoding algorithm and an error controler to compute a BER.

The basic architecture for demodulating two disctinct signals would be to duplicate this processing. Unfortunately, my hardware cannot contain so much computing.

One of the idea I had was to put in common some DSP operations. The downconversion uses a lot of ressources so I decided to mutualize the process for both channels.
The downconversion uses a CIC filter, it filters and decimates the signal.

The new signal path would be: DDS+MIXER (channel 1 brought to baseband, channel 2 centered around 50kHz), CIC + compensation filter (CIC response is not flat in passband), second dds + mixer to bring channel 2 to baseband.

here is a scheme of the processing:

The "tricky" thing is that on my second DDS I multiply both I and Q by cosine instead of cosine / sine (in that case, it does not work).
When doing this, I guessed that maybe I would lose something.

When measuring the BER, I can easily see that I have a 3dB loss between my two channels. My question: how can I "mathematically" explain this?

Is it because I am mixing two cosine instead of one cos and one sine; is half of the information lost?
Am I losing 3 dB on the SNR?

Edit: Here are the BER curves

Yes that is correct, you are only getting half the signal: with cosine on both inputs half the signal moves to baseband and the other half moves to twice the frequency offset. I suspect the reason it is not working in the sine cosine case is the sign of that is such that you translated the entire signal upward to twice the frequency offset instead of to baseband. If that is the case you could invert the sine input and it should work without the 3 dB loss.


How to re-program/modify sound frequencies without DSP processor? I'm developing a mobile application that controls Bluetooth sound-speaker functions, which includes an equalizer to modify the output audio frequencies.

The audio processing module on the speaker only has an Audio Codec processor (but no Digital Signal Processor, DSP). Is it possible to still implement a program that modifies the audio frequencies (but only relying on the DSP of the host device i.e Apple, Android device), and then transmits the modified frequencies to the speaker (Audio Codec Processor)? Thus, can you bypass the DSP?

Beginner

Is it possible to still implement a program that modifies the audio frequencies (but only relying on the DSP of the host device i.e Apple, Android device), and then transmits the modified frequencies to the speaker?

I can give you an absolute, positive, definite maybe.

You absolutely positively can do digital signal processing in your cell phone's CPU -- digital signal processing is just math, and can be done in any general-purpose processor.  The only thing special about a digital signal processor is that it's heavily optimized for the sort of computation that's needed for most DSP work.

The presence of dedicated DSP chips on the market is fading these days, because general purpose processors are getting cheaper and faster in general, and because general purpose processors are getting more and more "DSP-ish" features.  So there's a good chance that a new device will have plenty of raw horsepower.

Where the "maybe" comes in is whether the phone's operating system supports doing real-time DSP on the processor, whether the processor has sufficient computational power do implement your desired algorithm, and how easy the phone's operating system makes it for an app to do the work.  The answer to the OS part of that "maybe" is a question for Stack Overflow or some other more "pure" software-engineering oriented site; the answer to the processor loading part of that "maybe" is for you to determine experimentally once you've figured out how to write an app that can diddle with the audio stream.


How to reconstruct original signal from sampled signal? My original signal as
f1=2;
f2=5;

fs=100;
Ts=1/fs;
t=0:Ts:1;

xt=cos(2*pi*f1*t)+cos(2*pi*f2*t);

figure
plot(t,xt)

as shown below figure.

and my sampled signal with fs=12 which is greater than two times of
maximum frequency, 2*f2=10;
fs=12;
Ts=1/fs;
tn=0:Ts:1;

xn=cos(2*pi*f1*tn)+cos(2*pi*f2*tn);

and then I want to reconstruct the original signal with sampled signal.
I tried it using the matlab function interp1 as below
xr=interp1(tn,xn,t,'spline');

I compared original signal and reconstructed signal, but
they are different as shown below figure.

How to reconstruct original signal?

To implement a digital-to-analog converter, all you need is an ideal low-pass filter to filter out the periodic frequency response that $\varOmega \geq \varOmega_s/2$.
$$
H(j\varOmega) =\left\{
\begin{aligned}
T,\ \ \ \  |\varOmega|<\varOmega_s/2 \\
0,\ \ \ \  |\varOmega|\geq\varOmega_s/2
\end{aligned}
\right.
$$
The impulse response of an ideal low-pass filter is a sinc function
$$
h(t) = \frac{\sin(\pi t/T)}{\pi t/T}
$$
Convolving the sampled signal with this impulse response gives the analog signal.
$$
\begin{aligned}
y_a(t) =& \sum_{n=-\infty}^{\infty} x_a(nt) h(t-nT)\\
=&\sum_{n=-\infty}^{\infty} x(n) \frac{\sin[\pi (t-nT)/T]}{\pi (t-nT)/T}
\end{aligned}
$$
f1=2;
f2=5;

fs=100;
Ts=1/fs;
t=0:Ts:1;
xt=cos(2*pi*f1*t)+cos(2*pi*f2*t);
figure
plot(t,xt)

fs=12;
Tn=1/fs;
tn=0:Tn:1;
xn=cos(2*pi*f1*tn)+cos(2*pi*f2*tn);
hold on
stem(tn, xn)

y = dac(xn, Tn, t);
plot(t,y)

function y = dac(x, T, t)
m = 0:length(x)-1;
y = zeros(1, length(t));
for ii = 1:length(t)
    h = sinc((t(ii)-m*T)/T);
    y(ii) = sum(x .* h);
end
end


How to evaluate Spectral Density of DC component vs. AC component during computing PSD? When computing the PSD, how to manipulate the DC component?

PSD integral is the total signal power, does it include the power of DC component?

In other words, if there is one signal with DC value is high than 0, when we want to compute its PSD, do we need remove DC component first, then use function like pwelch() to compute PSD? If I want to evaluate the partition between spectral density of DC component and AC component, what should I do?

How to evaluate Spectral Density of DC component vs. AC component during computing PSD?
My final target is : Computing the partition between spectral density of DC component and AC component in one signal.

PSD is the spectral power density and also has the DC component, but in practice finite length of window inevitably produce leaking (because of side-lobes), so to prevent masking of small amplitudes near zero frequency we have to remove DC component.

Functions like pwelch don't remove DC component but it's better to remove it manually and after PSD calculations find it manually.


Precursor ISI - causality violation? Today, in my lecture on intersymbol interference (ISI), there was a discussion on postcursor ISI, and precursor ISI. Postcursor ISI is caused by the past bits, whereas, precursor ISI is caused by the future bits on the present bit. I could not understand how precursor ISI can occur. I was told that it is because of the delay of the channel. What does it mean?
Isn't there a causality violation? What am I missing?

Can someone please explain? Thanks.

Every symbol in the channel will be affected by ISI from previous symbols and next symbols. To understand this you must notice that there's a propagation delay from transmitter(TX) to receiver(RX). That means that it is possible for several symbols to be present in the channel at the same time before the first of them actually arrives at the receiver. When a symbol actually arrives at the receiver it will already have ISI from the symbols that were before and after it and shared that propagation time. Only the first symbol of the transmission (if the channel is free of previous symbols from perhaps some other packet, etc) will not have ISI from previous symbols. If we assume that the propagation time between TX and RX lasts, for instance what it takes to the TX to send 5 symbols then it is correct to assume that the 5 of them will interfere each other in more or less proportion because the are all present during the same time window in the channel. At the RX the first symbol will exhibit ISI from others which it hasn't received yet but were present in the channel.

You shouldn't analyze this as Oh, I'm receiving ISI from the future! Because the time t = 0 is not when the receiver received the first symbol but when the transmitter send it which gives you the whole propagation time to be causal.


Difference between Generalized Hough Transform and Cross-correlation Feature Matching I've been taking the Udacity Computer Vision Course, and one thing that has confused me so far is how the Generalized Hough Transform is any different from feature matching on an edge image using cross-correlation.

As I understand, the Generalized Hough Transform requires you to know the size, shape, and orientation of the object you are looking for, and using the edges of that object, parameterizes that object based on an R-table. I know it can be generalized to size and orientation, but that it is also costly to do. For now, let's say you choose not to generalize.

Given these constraints, couldn't you just do cross-correlation feature matching, using the same edge image the Hough Transform would be trained by?

What's the difference? Is there any reason Hough would perform better in this case, without generalizing to size and orientation?

There is a big difference: The Hough Transform maps the input space to a parameter space, where the search takes place. This way, the run-time of the algorithm is independent of the degree of the spatial search space. Correlation based methods are rather more brute force in that sense as they search explicitly for all transformations. Of course, there are ways to make that search more efficiently.

The Hough transform can easily be made rotation invariant at almost no cost. This generalizes to other types of linear transformations as long as the accumulator space doesn't over-grow as Hough tends to trade-off computational efficiency to storage.

Cross-correlation type methods are less memory intensive and are easier to extend to arbitrary transformations. Moreover, they allow for explicit control of the distance metric, giving rise to more creativity. Yet, they generally tend to be slow, unless implemented carefully. This is because they require a search over the whole input space for all (or part of) transformation space. Due to this this comparison paper, a well-implemented correlation approach with an efficient distance metric can outperform a Hough based one. Thus, we can say that many approximations to correlation methods can be made but one should make sure that the performance is not degraded.


Questions about Active Noise Cancellation Being new to ANC, this is my understanding of how ANC works:

Ambient noise (n) is entering a geographically small zone. To suppress the noise, an interfering (corrective) sound is injected acoustically into the zone by a speaker. The interfering sound from the speaker is adaptively generated so that it - in theory - ends up being a 180 degree phase-shifted version of the noise ....When this phase-shifted signal is acoustically mixed with the ambient noise, the noise is ideally cancelled out.

So based on the above, I would need a microphone to record sound as it presents itself in the "zone" and then let that drive an adaptive process which calculates a digital speaker signal to send to the speaker. The energy of the recorded sound would be the indicator of how well the adaptive process is performing (in terms of cancelling out the noise).

My question is : what type of algorithm(s) are used for this kind of DSP scheme?

Is one microphone enough?

What if the speaker has to output more than just the noise cancelling signal? For example: A user is wearing headphones with ANC. The user is listening to music but wants the "outside" noise to be cancelled out. I can't imagine that it would work unless the headphone has an "inside" microphone and an "outside" microphone.

Does ANC work in dimensionally small setups? For example:

The speaker is connected to a tube through which the corrective sound travels. The open end of the tube is placed in a small zone where a microphone is also located. An algorithm receives audio from the microphone, calculates the speaker signal and sends the speaker signal to the speaker. Will that work? Or are there some fundamental requirements/assumptions which have to be met for the ANC to work?

I would like to hear from people who wish to share their experience and knowledge acquired from working with ANC. Any advice is appreciated.

Thank you.

That's a pretty broad question. Let's start with a noise cancelling head set, that's about the easiest device.

A ANC headset has an internal microphone, that's placed as close to the ear canal as possible. It than runs a regular control loop (http://en.wikipedia.org/wiki/Control_theory): It calculates a speaker signal that tries to keep the signal at the microphone at zero. If you want to play music at the same time, you simply make the "desired signal" to be the music and not zero.

However, there are all sorts of implementation problems

Designing the control loop is difficult in terms of stability, cancelling performance, spectral fidelity, etc.
Manufacturing tolerances
Control loops need to be fast, so they work better at low frequencies. Most good headsets use a well tuned combination of active and passive cancellation. Passive works a lot better at high frequencies
While the headset cavity is small, it still has a fair amount of acoustic resonance and wave propagation, so you need to correct for the sound at the microphone vs. the sound at the ear drum
Many newer headsets use outside microphones as well, so it's a combination of feed-forward and feedback cancellation

Once the spaces get larger than a headset, things get a lot more complicated. Now you need to take into account the spatially distributed sound field, the environment, all kinds of propagation modes.


Audio Processing Glossary I'm new to audio processing and I hope here be a good place to ask about the definition of terms in this context.

My question is what is difference between sound, audio, voice, speech and tone, from both technical and mathematical point of view?

My own interpretation:

Sound: a mechanical wave that propagates through the air or water.

Audio: sound in the 20 Hz to 20 kHz range; in other words, sound that is (at least in theory) audible to a large number of humans.

Voice: sound produced by the human vocal tract.

Speech: intelligible voice (i.e. not grunts or screams)

Tone: signal dominated by a single frequency, even if it has harmonics.

Music: Harmonious, pleasant audio produced by man-made musical instruments. Often accompanied by song, which is audio just like music but produced by a human singer's vocal tract.

Musical tone: a tone produced by a musical instrument or singer.

Musical note: one of a particular, finite set of musical tones that, when played in combination, produce a pleasant auditory experience.


How does auto focus in cameras work? I am just wondering how does camera auto focus work. What kind of processing does it do on images in order to adjust the focus? Is there freely available code that detects the focus level of an image or video?

This answer actually provides free code to perform the focus detection.
https://stackoverflow.com/a/3795151/600736


Visualizing the measure of each FFT frame as an image in an audio signal i wrote a function in MATLAB which will divide a
audio waveform (a speech) into smaller windows with overlap .The output will be
a matrix X, each column of which contains the windows which i export.
Then i calculate the FFT of each column and i visualize the measure of each FFT frame as an image, using the "imagesc" command. I saw a shape that is symmetric
to the horizontal axis. Why i have this symmetry?.

Very simple:

The Fourier transform of every real-valued signal is (hermitian) symmetric.

You can thus, for your real-valued signals, simply not visualize half of the FFT bins.


Does the audio signal in the time domain change much after a RMS from one person to another? I would like to train a deep neural network to perform lip syncing given an audio input like in this article.

I want to train my neural network only on Obama speeches like done in the article and want to apply it to other people voices.
Does the structure of the audio signal in the time domain (after a RMS normalization) change much from one person to another if they say the same thing ?

In short, yes.

One thing that you'll find is that the human ear is actually relatively insensitive to phase, so you can take an audio sample, apply an allpass filter that doesn't affect the signal's magnitude spectrum but distorts its phase, and it won't make it sound much different to you. However, the time-domain waveform could look completely different. I'm not sure what techniques the researchers used in the link you gave, but I would guess that it isn't likely to be based upon the time-domain structure of the audio waveforms.


Time-domain LPF not showing expected behavior I am trying to implement a simple first-order Butterworth Low-Pass filter in Python. I have some code that makes use of scipy.signal.butter and scipy.signal.filtfilt. It works fine, but I wanted to make sure I could build the filter myself as well.

I start with a bilinear transform $s = \frac{2}{T_{s}}\frac{z-1}{z+1}$ on my first order LPF $\frac{\omega_{c}}{\omega_{c}+s}$. Making that substitution and converting from the $z$ domain to the time domain gives me the equation that I'm trying to reproduce:

$y_{n} = \frac{1}{\omega_{c}T_{s}+2}[\omega_{c}T_{s}(x_{n}+x_{n-1})+(2-\omega_{c}T_{s})y_{n-1}]$

I'm confident that this expression is correct, but my attempt at coding it doesn't seem to be working.

EDIT: I changed my implementation so what's actually happening is more clear. As a side effect, I'm now seeing a different type of strange behavior than I was previously. The new code is below:

T = 2                                                  # time length over which samples are collected (s)
f0 = 5                                                 # fundamental frequency (Hz)
fs = 1000                                              # sampling frequency
Ts = 1/fs                                              # sampling period
fc = 10                                                # filter cutoff frequency
wc = 2*pi*fc                                           # cutoff frequency in radians

# Create wave input to be sent to filter

t = np.linspace(0,T,T*fs+1)
x_storage = [10 + cos(2*pi*f0*x) for x in t]
plt.plot(t,x_storage)

# Create a place to store output data

y_storage = np.zeros(len(t))

# Initialize input and output data

y = 0                                                  # y_(n)
yprev = 0                                              # y_(n-1)
x = 0                                                  # x_(n)
xprev = 0                                              # x_(n-1)

# Begin calculation

for idx in range(len(t)):
    # One timestep of the calculation
    xprev = x                                          # move the old x_(n) value to x_(n-1)
    x = x_storage[idx]                                 # bring in the current datapoint to be the new x_(n)
    yprev = y                                          # move the old y_(n) value to y_(n-1)
    y = (wc*Ts*(x-xprev)+(2-wc*Ts)*yprev)/(2+wc*Ts)    # calculate the new y_(n) according to the filter equation
    y_storage[idx] = y                                 # store the y_(n) value for later plotting

 plt.figure()
 plt.plot(t,y_storage)

My input is a wave with a DC offset and a frequency component below the cutoff threshold, so I would expect the output wave to look more-or-less the same. Instead, I see the DC bias vanish and the signal itself get hugely attenuated (see below). Any thoughts would be much appreciated!

Your math is correct but your code doesn't match it:

y = (wc*Ts*(x-xprev)+(2-wc*Ts)*yprev)/(2+wc*Ts)

It should be $x+xprev$ not $x-xprev$. The minus sign makes it a highpass and hence the attenuation.


Digital Audio Amplifiers Can anyone tell me why a filter that blocks DC signal in an Amplifier would have a cutoff frequency of 4 Hz . Why would someone choose a frequency as low as 4Hz to filter out DC content before processing the audio signal further ?

In order to block the DC component in a signal, the high-pass filter can be used.

In the figure above, a basic high-pass filter which is known as the passive RC high-pass filter can be seen. The working principle of it depends heavily on the capacitor’s behaviour at varying frequencies as the formula depicts:
$$X_{C}(j\omega) = \frac{1}{j\omega{C}} \Omega$$
The formula is named as the capacitive reactance which is actually a capacitor’s resistance (not the internal resistance) by which an alternating current is attenuated while it passes through it.
So, let’s analyse the given circuit. If the Kirchhoff’s voltage law is applied on that circuit;
$$-x(t) + {\frac{1}{C}\int_{t_{0}}^{t} i(t) \, dt} + V_{C}(t_{0}) + i(t)R = 0$$
$$-x(t) + {\frac{1}{C}\int_{0}^{t} i(t) \, dt} + V_{C}(0) + i(t)R = 0$$
Let’s assume that the capacitor’s inital voltage $V_{C}(0) = 0V$.
$$-x(t) + {\frac{1}{C}\int_{0}^{t} i(t) \, dt} + i(t)R = 0$$
$${\frac{1}{C}\int_{0}^{t} i(t) \, dt} + i(t)R = x(t)$$
If the last equation is differentiated with respect to time and the relation $y(t) = i(t)R$ is used according to the circuit;
$${\frac{1}{C}i(t)} + R\frac{di(t)}{dt} = \frac{dx(t)}{dt}$$
$${\frac{1}{RC}y(t)} + \frac{dy(t)}{dt} = \frac{dx(t)}{dt}$$
By rearranging the last equation, we obtain;
$${\frac{dy(t)}{dt} + \frac{1}{RC}y(t)}  = \frac{dx(t)}{dt}$$
Interpreting what is happening in a high-pass filter, low-pass filter, etc. by looking at their differential equations can be demanding. Since these circuits are used to make changes on signals with regards to changes in frequencies, analysing these circuits in the frequency domain can be quite useful. In order to do this, the Fourier transforms of the both sides of the last differential equation can be taken.
$${j\omega}Y(j\omega) + {\frac{1}{RC}}Y(j\omega) = {j\omega}X(j\omega)$$
$${Y(j\omega)}(j\omega + \frac{1}{RC}) = {j\omega}X(j\omega)$$
Consequently, the transfer function $H(j\omega)$ becomes;
$$H(j\omega) = \frac{Y(j\omega)}{X(j\omega)} = \frac{j\omega}{j\omega + \frac{1}{RC}}$$
As it can be seen from the formula, basically, the transfer function results in a complex number. In order to obtain its magnitude, the absolute value of it can be calculated.
$$|H(j\omega)| = |\frac{Y(j\omega)}{X(j\omega)}| = |\frac{j\omega}{j\omega + \frac{1}{RC}}|$$
$$|H(j\omega)| = \frac{|j\omega|}{|j\omega + \frac{1}{RC}|}$$
$$|H(j\omega)| = \frac{\omega} {\sqrt{ {\omega^{2}} + \frac{1}{R^{2} C^{2}} }}$$
This equation relates the input voltage of the passive RC high-pass filter to the output voltage by the input voltage’s frequency. Actually, this equation explains the main working principle of the passive RC high-pass filter by itself.
In order to make the transfer function’s behaviour i.e. the high-pass filter’s behaviour over frequency perceptible, let’s create an example circuit with defined values of $R$ and $C$ and illustrate the magnitude graph of the transfer function. I’ve chosen the values arbitrarily and assumed that $R = 20 k\Omega$ and $C = 1 \mu{F}$.
The MATLAB code for drawing the magnitude graph of the transfer function is given below:
f = 0: 0.01: 20;  % Frequency range in Hz

R = 20000;  % Resistance value in Ohms
C = 0.000001;  % Capacitance value in Farads

omega = 2 * pi * f;  % Angular frequency in rad/s

mtf = omega ./ sqrt((1 / power(R * C, 2)) + power(omega, 2));

plot(f, mtf, "r", "LineWidth", 2);
grid on;
xlabel("Frequency (Hz)");
ylabel("|H(f)|");
title("Frequency Response of the Passive RC High-Pass Filter");

I’ve plotted the graph with respect to the frequency as most of the time, we are dealing with classic frequency rather than angular frequency.
The resultant plot is given below:

As it can be seen from the plot above, the corner/cut-off frequency of the example high-pass filter is around 8 Hz which can be calculated by the formula;
$$f_{C} = \frac{1}{2\pi{RC}} Hz$$
Corner/cut-off frequency is the frequency at which the output signal’s power which is transferred to the load (in our case, it is the resistor $R$) is half of that of the input signal. Alternatively, it is the point at which the output voltage’s magnitude is approximately 0.707 times the input voltages’s magnitude.
So, for our passive RC high-pass filter; if the input signal’s frequency is chosen lower than the corner/cut-off frequency of the filter, the input signal will be attenuated dramatically.
Since a DC signal has the frequency of $f_{DC} = 0 Hz$ ideally, it can be removed from the input signal.
However, I don’t have any idea about why the manufacturer has chosen the cut-off frequency as 4 Hz but choosing it lower than the minimum frequency value of the human hearing range which is 20 Hz is obvious as we don’t want to miss a frequency component while we are listening to music, conversations, etc. via a loudspeaker.
This specific cut-off frequency rating could be about the designer’s state of mind or the circuitry’s needs.


Linear Panorama Transform I am trying to create a linear panorama (a series of images taken from the same plane but not the same location and not rotated) but was not sure which transform is appropriate for this task.  I know a homography is used to stitch images together from a rotated point.  I think a transform is necessary but I am not even sure of that.

No transform is required for that kind of thing except making sure that the rays of light hit the lens head on.

This is only possible when the subject is far away enough from the camera and for a narrow range of its field of view.

In other words, the end result image we are talking about here is composed of vertical projections of 3D points right on the back of the camera. As if you had a one pixel lens which you scanned across the full apperture of your scene to sample exactly what is right at the focus point of that lens. Head on.

But what really happens is that for every picture you take, the light goes through a single lens. The only pixel that is vertically sampled (or projected, from the world to the sensor) in this case is the one at the centre of that lens. Everything else is imaged at an angle that is proportional to how far away from the center of the lens the sampling pixel is located.

Take the picture of a very tall multistorey building for example. At eye level, straight ahead, the windows appear rectangular. But as you raise your gaze upwards, you "sample" the windows higher up at an increasing angle to the point that some of the windows might even start to appear square (their perceived height, approximately equal to their perceived width).

If you do not have a 3D mesh of the scene, there is no transform that will correct for this change in perspective.

So, what you do is to "fake it".

You crop a patch from the acquired images that given the size of the subject, its distance from the camera and the focal length of the lens, it appears to be imaging the subject head on. And you take lots of images.

No further reprojection is required except of course to correct for the relative camera poses between takes as the camera is not on a leveled rail but that's scaling and rotation only.

Depending on the application this is sometimes called the Pushbroom camera model.

Provided that you have taken care during acquisition, Hugin can do that type of stitching.

Hope this helps.


$n\textrm{ Hz}$ waveform sampled at $m\textrm{ Hz}$ per second Here is an example of plotting a square wave given in SciPy Documentation

A $5\textrm{ Hz}$ waveform sampled at $500\textrm{ Hz}$ for 1 second:

from scipy import signal
import matplotlib.pyplot as plt
t = np.linspace(0, 1, 500, endpoint=False)
plt.plot(t, signal.square(2 * np.pi * 5 * t))
plt.ylim(-2, 2)**strong text**

The documentation page states that

  The square wave has a period 2*pi, has value +1 from 0 to 2*pi*duty
  and -1 from 2*pi*duty to 2*pi. duty must be in the interval [0,1].

Why in the example 5*t is multipled by 2*np.pi?
Doesn't the signal.sqrt() function
take care of multiplying the frequency 5*t by 2*np.pi?
What is the meaning of sampled at $500\textrm{ Hz}$ for second? I understand it is a $5\textrm{ Hz}$ waveform since 2*np.pi*5*t takes care of creating $5\textrm{ Hz}$ wave. But I do not understand

  sampled at 500 Hz for 1 second

The basic square wave function signal.square() has a period of 2*pi. Multiplying each time by 2*pi*f is a way of stretching the wave along the time axis to get a wave with a period of 1/f.
To understand what the sample rate refers to, look at the first line:

t=np.linspace(0,1,500,endpoint=False)

If you look at the numpy documentation for this function, you'll see this creates a list of 500 evenly spaced numbers between 0 and 1. This is an array of times. When you then call the signal.square() function with this array, it returns an array of points on the square wave for each of those times. Thus you have sampled the square wave at 500Hz for 1 second.


Why does block LMS have the same performance as LMS? The block LMS and conventional LMS have the same convergence rate and the same misadjustment. I am having trouble wrapping my head around this. The block LMS uses a more accurate estimate of the gradient vector at each iteration. Conceptually, why does a better gradient estimate have no benefit on the descent?

here is a short and little (a.k.a. "consise but terse") derivation of the LMS and normalized LMS adaptive filter.

once the LMS has converged on a reasonably stable equilibrium for $h_n[k]$, they won't move around that much.  then it doesn't matter so much how long the block is.  and the only difference between block-LMS and the plain-old ordinary LMS is the block size.  (the block size for the latter is 1.)


What is the difference between Spectrogram and Scalogram? Could someone please explain the difference between the two?

Spectrogram:

  A visual representation of the spectrum of a sound changing through time.

Scalogram:

  (signal processing) A visual representation of a wavelet transform, having axes for time, scale, and coefficient

Both convert a 1D signal into a 2D complex map, aiming at unveiling transient features of the data. Both use a functional microscope: the signal $s(t)$ is analyzed through variations of a single template function $\rho(t)$ with two parameters $p$ and $q$.

So the signal is analyzed through the form:

$$R_s(p,q) = \int s(t) \rho^*_{p,q}(t) d\cdot$$

which is 2D (the differential increment $d\cdot$ can be specific).

With the spectrogram, $\rho$ is a fixed (non-zero sum) window, and $\rho_{p,q}$ is the same window, shifted by $p$, and  modulated by a frequency proportional to $q$.
With the scalogram, $\rho$ is a fixed zero-mean function (a wavelet), and $\rho_{p,q}$ is the same wavelet, shifted by $p$, and  dilated by a factor proportional to $q$.

Why they have those properties is another story.


Converting magnitude to dBs I have seen people using 10*log10 and/or 20*log10 while converting magnitude to dBs. What are the differences and which one is valid?

They are both valid:

$10 \cdot \log_{10}$ is for power ratios
$20 \cdot \log_{10}$ is for voltage/pressure/etc ratios (these are called "field quantities" or "root-power quantities")

See Field, power, and root-power quantities and Decibel definition

Digital values are just unitless numbers, so you need to think about what they ultimately represent.  Typically they're root-power quantities (like audio pressure waveforms), but could conceivably represent power waveforms.


Why ups and downs of amplitude? Ladies, Gentlemen,
Please let me again employee you with my problems. In following program there is only one sinusoid with frequency 256 Hz amplitude 100,000 and phase 0.67 pi rad. Testing frequencies grid is of 1 Hz step, from 0 Hz to 1000 Hz. Output has maximum amplitude 99,221 and phase 0.670561 pi rad at 256 Hz. I regard it okay given errors in accuracy.

But there are ups and downs in pseudo-frequencies amplitude while my understanding about DFT is that they should going upward till 256 Hz and then downward. Why these ups and downs? I remain with regards.

// Grid is created by variables lowTestFreq and testFreqStep
// both of double type, and integer type M.
// In comments are units of variables.

#include <iostream>
#include "math.h"

int main()
{
    //Unchanged variables.
    double pi = 3.141592653589793; //pure number.
    int n = 0; // sample.
    double xn = 0; //pure number.
    double Re = 0; //pure number.
    double Im = 0; //pure number.
    double xn_cos = 0; //pure number.
    double xn_sin = 0; //pure number.
    double testAmpl = 0; // pure number.
    double testPhase = 0; // rad and then pirad.

    //Sampling frequency and samples of segment.
    int samplFreq = 50000; // sample/second.
    int N = 3000; //sample.

    // Existing sinusoid.
    double ampl_1 = 100000; //pure number.
    double existFreq_1 = 256; // Hz.
    double phase_1 = 0.67 * pi; // rad.
    double unitAngle_1 = (existFreq_1 / samplFreq) * (2 * pi); // rad/sample.

    // Testing grid.
    double testFreq = 0; // Hz.
    double testUnitAngle = 0; // rad.
    double testAngle = 0; //
    double lowTestFreq = 0; //Hz.
    double testFreqStep = 1; // Hz.
    int M = 1000; // pure number.
    int m = 0; // pure number.

    for (m = 0; m <= M; ++m)
    {
        testFreq = lowTestFreq + (m * testFreqStep); // Hz.
        testUnitAngle = (testFreq / samplFreq) * (2 * pi); // rad/sample.

        Re = 0; Im = 0;

        for (n = 1; n <= N; ++n)
        {
            xn = ampl_1 * sin(n * unitAngle_1 + phase_1);

            testAngle = n * testUnitAngle;
            xn_cos = xn *  cos(testAngle);
            xn_sin = xn * -sin(testAngle);

            Re += xn_cos;
            Im += xn_sin;
        }
        std::cout << "Test Freq " << testFreq; // Hz.
        testAmpl = sqrt(Re * Re + Im * Im) / N * 2; //pure number.
        std::cout << "       Test Freq's Ampl: " << testAmpl;
        testPhase = atan2 (Im, Re); // rad.
        testPhase = testPhase / pi; // pirad.
        testPhase = testPhase + 0.5; // pirad. Corrected phase.
        std::cout << ",    and its phase: " << testPhase
            << " pirad." << std::endl;
    }
    return 0;
}

A DFT includes an implicit finite length window on any function (which has an infinite domain in the case of a sinusoid).  The effect of multiplication by that window in the time domain on the DFT result is a convolution with the transform of that window.  In the case of a rectangular window, that looks like something very close to a Sinc function (except near DC or Fs/2), or sin(w)/w (which has lots of "ups and downs".)


Reducing number of taps in FIR for echo cancellation I'm trying to implement the NLMS algorithm on an FPGA. The problem is that I'm only allowed to use a maximum number of 2048 coefficients for my FIR filter, my sampling rate is 16KHz, however, the room impulse response sampled at 16KHz can be 10 times longer than 2048, is there any way to decrease the number of taps to remove long echos with only 2048 taps? I'm new to FIR subject, I know that my FIR filter IP allows decimation, interpolation, and can have different sampling rate than the input, can I exploit this to decrease the number of taps? I can also compress room impulse response and round many samples that are close to zero and make them zero but can this help to reduce the number of coefficients?

The only data on the environment where your setup is operating is given in the phrase "the room impulse response sampled at 16KHz can be 10 times longer than 2048". These 20480 samples cover a time interval of 1280 ms. I wonder what the room it is, because in the article "Adaptive algorithms for sparse echo cancellation" I read: "The length of the acoustic echo response in a typical teleconferencing room is in the region of 100 to 400 ms". The good news for you is in the phrase continuation: " and hence adaptive filters employing 1024 taps or more are typically required in order to achieve adequate levels of echo cancellation". Therefore, using the techniques described in the ref., with your 2048-tap filter you will be able to cancel echo in the region of up to 800 ms, not far from claimed 1280 ms.
I recommend you first to make sure that your estimation (and understanding) of the room impulse response parameter is correct. If the task specs are more optimistic, you have the ready solution among the adaptive algorithms of the article "Adaptive algorithms for sparse echo cancellation". If the task specs are tougher, your best bet is to use further developments of adaptive filter algorithms, as a nonparametric variable step-size PNLMS algorithm (NPVSS-NLMS) and a sparseness-controlled affine projection algorithm (PAPA).
Being a novice in the subject, you may benefit from learning of a first-entry approach to a sparse array-based room transfer function estimation for echo cancellation.


How does the $\mathcal Z$-transform's "region of convergence" work? I'm a novice in DSP and I have few doubts regarding the $\mathcal Z$-transform and its region of convergence (ROC).

I know what a $\mathcal Z$-transform is. But I'm having trouble with understanding the ROC. First of all I have some confusion with $X(z)$ and $x(z)$. I easily get caught by exchanging these terms. I know the ROC defines the region of where the $\mathcal Z$-transform exists. From the web and my books states that:

  If $x[n]$ is a finite-duration sequence, then the ROC is the entire $z$-plane, except possibly $z = 0$ or $\lvert z\rvert = \infty$. A finite-duration sequence is a sequence that is nonzero in a finite interval $n_1 \le n \le n_2$

And later it says:

  When $n_2 > 0 $ there will be a $z^{-1}$ term and thus the ROC will not include $z=0$. When $n_1 < 0$ then the sum will be infinite and thus the ROC will not include $\lvert z\rvert=\infty$.

This is where I get stuck!. What they try to say in the above line "When $n_2 > 0$ there will be a $z^{-1}$ term and thus the ROC will not include $z=0$" What do they mean by $z=0$? Are they substituting $z$ as $0$, if so in which equation?

How do we calculate the region of convergence for an infinite sequence?

To be completely honest, I thought the theory behind the Z-transform was kinda opaque in college too. In hindsight, taking a course in complex analysis would have made it clearer. And I too dislike the notational conventions that seem to be used for this stuff. Strictly speaking, the usual convention here is that

$x[n]$ denotes a discrete-time sequence

$n\in \mathbb{Z}$
the brackets denote a discrete argument

$X(z)$ denotes a continuous-valued transformed function

$z\in\mathbb{C}$ (is a complex number)
the parentheses denote a function accepting a continuous-valued parameter
the capital $X$ denotes a transformed version of some other function/sequence $x$ (a similar notation is used for Fourier transforms: $F(j\omega)\leftrightarrow f(t)$

  What do they mean by z=0? Are they substituting z as 0, if so in which equation?

They mean, just plug $z=0$ into your usual definition of the Z-transform.

$X(z) = \sum_{n=\infty}^{\infty} x[n] z^{-n}$

Generally (more precisely, when $x[n] \ne 0$ for some $n \ne 0$), this sum will diverge (to infinity) for some complex $z$. For instance, let $x[0] = 1, x[1] = 1$, and $x[n]=0$ for $n<0$ and $n>1$. Then $X(z) = 1 + z^{-1}$. The ROC does not includes $z=0$, for $\lim_{z\rightarrow 0} X(z)=\infty$

When your text says "When $n_2>0$ there will be a $z^{−1}$ term and thus the ROC will not include $z=0$", what they mean by that is, when $x[n]$ is nonzero for some $n>0$, it is unavoidable for the z-transform to include a $z^{-n}$ term, which diverges to infinity at $z=0$. That's all.

  How do we calculate the region of convergence for an infinite sequence?

Lots of math. Ha!

srsly, the way this is done is to obtain an algebraic formulation for the sequence in question, plug it into the Z-transform definition, and use the tools available from the analysis of geometric series (and complex power series) to determine where this Z-transform converges/diverges. In practice, determining whether $|z|=1$ converges is the most important question to answer, because that determines stability, and whether or not you can obtain a frequency response from the system, etc. But causality might matter too, depending on what you're doing.


Autocorrelation Jitter I'm doing autocorrelation (via fft) on pure, synthetic generated sinusoids and getting unexpected jitter in the results.  The jitter depend on where in the sine the signal buffer is started.  I get the correct results if I start the buffer at 0 and 180 degrees of a sine wave, and get the most error at 90 and 270 degrees.  I am padding with zeros after the signal buffer up to the next power of 2 above 2x-1 the signal buffer size.

I've been using scipy fftconvolve, as well as doing the auto-correlation with fft, multiply the results by it's complex conjugate, then inverse fft.   I get the same results both ways.   This makes me think I am not preparing my signal buffer properly or analyzing the autocorrelation function results correctly.

Below are two examples of autocorrelation of a 440hz signal.  The signal buffer the autocorrelator is analyzing is enclosed by the two vertical lines, approximately two and a half periods and 256 samples at 44.1khz sample rate.

In the first example, the buffer is filled with the sine starting at 0 degrees.  The result is 440.7 hz.

In the second example, below, the buffer is filled with the sine starting at 90 degrees.  The result is 453hz, a considerable error.

Any ideas on why the error is so large in the second example?

The (periodic) autocorrelation function of a sinusoid $A\cos(2\pi f_0 t + \theta)$
is

$$\begin{align}
R(\tau) &= \int_0^{f_0^{-1}} A\cos(2\pi f_0 t + \theta)
A\cos(2\pi f_0 (t +\tau)+ \theta)\,\mathrm dt\\
&= \frac{A^2}{2}\int_0^{f_0^{-1}} \cos(2\pi f_0\tau)+\cos(2\pi f_0 (2t +\tau)+ \theta)\,\mathrm dt\\
&= \frac{A^2}{2f_0}\cos(2\pi f_0\tau)
\end{align}$$
which has a peak at $\tau = 0$. Notice in
particular that there is no dependence on $\theta$.   More generally,
if the integration is carried out over an interval $[0,T]$
of arbitrary length, we get
$$\begin{align}
\hat{R}(\tau) &= \int_0^{T} A\cos(2\pi f_0 t + \theta)
A\cos(2\pi f_0 (t +\tau)+ \theta)\,\mathrm dt\\
&= \frac{A^2}{2}\int_0^{T} \cos(2\pi f_0\tau)+\cos(2\pi f_0 (2t +\tau)+ \theta)\,\mathrm dt\\
&= \frac{A^2T}{2}\cos(2\pi f_0\tau) + \left.\frac{A^2}{2}\frac{\sin(2\pi f_0 (2t +\tau)+ \theta)}{4\pi f_0}\right|_0^T\\
&= \frac{A^2T}{2}\cos(2\pi f_0\tau) + \frac{A^2}{2}\frac{\sin(2\pi f_0 (2T +\tau)+ \theta)-\sin(2\pi f_0\tau + \theta)}{4\pi f_0}\\
&= \frac{A^2T}{2} \left[\cos(2\pi f_0\tau) + \frac{\sin(2\pi f_0 (2T +\tau)+ \theta)-\sin(2\pi f_0\tau + \theta)}{4\pi f_0T}\right]
\end{align}$$
which does depend on $\theta$ as well as $T$. Notice that
if $T$ is an integer multiple $kf_0^{-1}$ of the fundamental
period $f_0^{-1}$, then the second term
in brackets disappears entirely and we get
$$\hat{R}(\tau) = kR(\tau).$$
On the other hand, since the numerator in that second term cannot
have value more than $2$, we get that if $T \gg f_0^{-1}$ but
is not an integer multiple of $f_0^{-1}$, the
second term is quite small in magnitude and so the autocorrelation
function that one deduces from $\hat{R}(\tau)$ depends only
very slightly on $\theta$.

With that as prologue, I point out to you that the FFT method of
computing the autocorrelation function of a signal gives the
periodic autocorrelation function of a segment of the
signal, that is, what one gets is effectively an approximation to
$\hat{R}_x(\tau)$ for a periodic signal $x(t)$ that is approximately
$2.5$ periods of a sinusoid repeated over and over again.  In particular,
this periodic repetition is not the same as the sinusoid continuing
smoothly and without break as one would want it to be for the above
results to apply.  Thus, it is not at all surprising that the
results from the
FFT method of computing the autocorrelation function depend more
heavily on the phase of the sinusoid. In approximately
$2.5$ periods, the $sin(\cdot)$ goes from $0$ at $t=0$ to nearly
$0$ again at the end of nearly $2.5$ periods. Thus, the periodic
extension is nearly continuous (though there is an abrupt
change of phase at the transition). On the other hand, with
a $90$-degree phase shift, $\cos(\cdot)$ goes from $1$ at
$t = 0$ to nearly $-1$ at the end of nearly $2.5$ periods,
and so the periodic extension is very discontinuous. This
shows up in the FFT results which don't get a good estimate of
the frequency in the latter case.

If you take the frequency of the sinusoid to be be such
that an integer number of cycles fit into the $N$ that
you are using for your FFT, these artifacts will disappear.


Detecting and isolating part of an image I'm working with cross-sectional CT scans of the leg using matlab's image processing toolbox.  Part of what I'm attempting to do is to automatically isolate the leg in any image, which requires removing any extraneous objects that are present in the image.  So far my code is able to take out most objects, which works for the majority of my images.  However, some of the images have an object that is flush with the leg, which my code isn't able to remove.

Here's an example of the problem: The first image shows the original image, and the second one shows the "filtered" image, after it has gotten rid of extraneous objects.  But as you can see, there is still something outside of the leg cross-section that I would like to get rid of.

I've thought of using something to detect ellipses, but unfortunately some of my images are not very good homogenous shapes (like the one below).

And here's the current code I have:

% Read input file
inputfile = dlmread(filename,';');

% Normalise input file (0 to 1) & show original image
I_original = mat2gray(inputfile,[min(inputfile(:)) max(inputfile(:))]);
figure; imshow(I_original)

% Detect leg
BW1 = im2bw(I_original,.2);
BW2 = imfill(BW1,'holes');

% Remove misc. objects around leg (i.e., stray objects <5000 pixels)
BW_final = bwareaopen(BW2,5000);

% Generate filtered image
I_filtered = I_original;
I_filtered(imcomplement(BW_final)) = NaN;
figure; imshow(I_filtered)

Any suggestions on how to fix this would be appreciated!  Thanks!

Since the Extraneous part is connected with leg section via a small thread like thing, Try converting it to binary Image, and then eroding the image to get rid of smaller objects. That should disconnect your both objects. Then choose the object with biggest area.

Helpful functions - im2bw, imerode, rogionprops.


2-D Wavelet output changing when using 'figure' command in Matlab I am having trouble understanding what is causing the problem to my Wavelet output when i am trying to display the image with the figure command.

I want to display both the Original image and its Wavelets decomposition,
but when I'm using this code the original image is not displaying in a separate figure.

Here is my code,

im = imread('cameraman.tif');
imshow(im),title('Original Image')

[cA1,cH1,cV1,cD1] = dwt2(im,'bior3.7');
A1 = upcoef2('a',cA1,'bior3.7',1);
H1 = upcoef2('h',cH1,'bior3.7',1);
V1 = upcoef2('v',cV1,'bior3.7',1);
D1 = upcoef2('d',cD1,'bior3.7',1);

% To display the results of the level 1 decomposition, type:

subplot(2,2,1); image(wcodemat(A1,192));
title('Approximation A1')

subplot(2,2,2); image(wcodemat(H1,192));
title('Horizontal Detail H1')

subplot(2,2,3); image(wcodemat(V1,192));
title('Vertical Detail V1')

subplot(2,2,4); image(wcodemat(D1,192));
title('Diagonal Detail D1')

and the output of this code is just one figure with four plots (Digonal, Horizontal, Vertical details and Approximations)

but when I'm using figure command to display the original image in a separate figure,

 im = imread('cameraman.tif');
imshow(im),title('Original Image')
figure

it messes up the wavelet output,

Why does it happens that my output is changed to blue and red color ?

The reasoning is simple. MATLAB defaults to a 'jetblue' colormap when you use imshow(im).

Typically you need to specify the colormap being used when you use imshow.

For this, add another command like:

$$colormap(gray(256));$$

to your current script.


How to remove frequency offset with RRC as pulse shaping filter I am working on a project with QPSK modulation format and a RRC filtering at the transmitter (TX) side. The channel introduces a constant dropper shift ($f_d$) to the received signal. To make the channel ISI free, the same RRC filter is used on the receiver (RX) side.
To remove this $f_d$, a practical ML frequency offset estimator is used [1].
This method is derived by assuming this $f_d$ is small enough. Therefore, after some mathematical manipulation, the decision statistics can be narrowed down to
$$\max(\lvert \text{fft}(c(k)^**y(k)\rvert),$$
where y(k) is the received r(k) filtered by a RRC matched to the RRC at TX, and c(k) is the transmitted data symbol (sample). To eliminate c(k) in the statistics, the $(c(k)^**y(k) )^4$ is used before fft operation is performed. The upsampling rate is 10. I couldn't see any issues in this process. Besides, I also tried to do $((c(k)*y(k))^4 )^{\frac 14}$ in order to eliminate the effect on the magnitude of the filtered received signal, but this gives me a huge DC component.
Finally, the algorithm works very well without RRC filtering at TX and RX. It also works very well without this $f_d$. So the problem should be something related this pulse shaping part.
Could anybody see any problems in this process?
[1]. U. Mengali et al: "Synchronization Techniques for Digital Receivers".

The RRC output only takes the exact constellation point values at exactly the right timing instant. That's what I meant with "coherent detector feeding a synchronously sampled matched filter" in an answer to a previous question of yours.
You'll have to deal with the inter-symbol interference with previous symbols. Which is no big deal – all symbols should have the same probability, so the average phase error you'd get is zero.
But, not so much for the fourth power: that's always a positive error, so yes, you get a biased phase correction term, and that means a misestimate of your frequency. Although I only mentioned noise, uncorrelated ISI looks a lot like nose, and thus it should have been among the reasons when I said

I'd advise against this whole approach!

in the answer to a question that Linda posted a couple of months ago (and seeing that this deleted question of four hours ago is practically a precursor to this question, I'll assume you and Linda basically share a team, or actually are two accounts for the same person. If the second account was made by mistake: Moderators can merge these for you, no problem! Just ask over in Meta.DSP).
You can now either:

Do the timing recovery before the frequency correction, and live with the significantly reduced performance that brings due to ISI or
Drop the fourth-power approach and go for a second-order PLL.


Received Signal Error vs. BER It is my understanding that the least squares algorithm (e.g., in equalization) minimizes the received signal error. However, minimizing the received signal error does not necessarily equate to minimizing the BER.

I was wondering why those two are not always equivalent? And how to properly gauge the performance of the least squares algorithm if BER is not the right metric?

EDIT: Simulation description

random bits --> frame bits --> NRZL --> modulator --> AWGN --> demodulator (integer multiple oversampled) --> LS EQ --> soft bits --> BER comparison

modulator: FSK and BPSK (tried both)
AWGN: Eb/No of 7dB - 20dB
frame structure: 16-bit training sequence followed by 496 random bits
Optionally add multipath delay as well (1-2 bauds).
Taps on LS EQ 3-11 bauds.

I will try to answer on your question with this: Least Square equalizer is minimizing the intersymbol interference of your signal (ISI), between the symbols, which does not translate directly to minimization of the bit error rate. When I am saying 'directly' that means that there is no exact linear dependency between the error reduction in the equalizer, and the bit error rate. However, they generally tend to go into the same direction.

The reason for that is the nature of equalizer, as the nonlinear estimator.

The chain of processing that you have shown consist of processing modules, and some of them are nonlinear, which interferes with the bit error rate  minimization. In other words, bit error rate minimization is not a criterion that is used by the equalizer.

Below is the mass of papers that are talking about it, and it is not a simple matter.

Please note the first one: "Approximate Minimum Bit-Error Rate Equalization for Binary Signaling" which actually argues exactly the same line of inquiry that you are advocating: "Although most linear and decision-feedback equalizers are designed to minimize a mean-squared error (MSE) performance metric, the equalizer that directly minimizes biterror rate (BER) may signiﬁcantly outperform the minimum MSE equalizer".

This means that the digital bit error based criterion would be a better choice for reducing the bit error rate than the one that is typically used.

http://www.ijesit.com/Volume%202/Issue%203/IJESIT201303_07.pdf

https://www.mathworks.com/help/comm/examples/ber-performance-of-different-equalizers.html?s_tid=gn_loc_drop

http://barry.ece.gatech.edu/pubs/conference/icc97.pdf


What algorithm does Google use for its "Search By Image" site? What is your best guess how Google Image Search works? I can upload a photo and can search for similar images. What algorithm does it use to identify similar images?

I don't know which algorithm Google uses. But, since you wanted a best guess, let me give some ideas on how a similar system could be constructed.

The whole field dealing with search-image-base-by-image is called Content Based Image Retrieval (CBIR). The idea is to, somehow, construct an image representation (not necessarily understandable by humans) that contains the information about image content.

Two basic approaches exist:

retrieval using low-level (local) features: color, texture, shape at specific parts of images (an image is a collection of descriptors of local features)
semantic approaches where an image is, in some way, represented as a collection of objects and their relations

Low-level local approach is very well researched. The best current approach extracts local features (there's a choice of feature extraction algorithm involved here) and uses their local descriptors (again, choice of descriptors) to compare the images.

In newer works, the local descriptors are clustered first and then clusters are treated as visual words -- the technique is then very similar to Google document search, but using visual words instead of letter-words.

You can think of visual words as equivalents to word roots in language: for example, the words: work, working, worked all belong to the same word root.

One of the drawbacks of these kinds of methods is that they usually under-perform on low-texture images.

I've already given and seen a lot of answers detailing these approaches, so I'll just provide links to those answers:

CBIR: 1, 2
feature extraction/description: 1, 2, 3, 4

Semantic approaches are typically based on hierarchical representations of the whole image. These approaches have not yet been perfected, especially for the general image types. There is some success in applying these kind of techniques to specific image domains.

As I am currently in the middle of research of these approaches, I can not make any conclusions. Now, that said, I explained a general idea behind these techniques in this answer.

Once again, shortly: the general idea is to represent an image with a tree-shaped structure, where leaves contain the image details and objects can be found in the nodes closer to the root of such trees. Then, somehow, you compare the sub-trees to identify the objects contained in different images.

Here are some references for different tree representations. I did not read all of them, and some of them use this kind of representations for segmentation instead of CBIR, but still, here they are:

binary partition trees and mention of min/max trees: P. Salembier, M.H.F Wilkinson: Connected Operators
binary partition trees: V. Vilaplana, F. Marques, P. Selembier: Binary Partition Trees for Object Detection
tree of shapes (component tree): P. Monasse, F. Guichard: Fast Computation of Contrast-Invariant Image Representation, C. Ballester, V. Castellis, P. Monasse: The Tree of Shapes of an Image
monotonic trees: Y. Song, A. Zhang: Analyzing scenery images by monotonic tree
edit: further digging shows that tree of shapes and monotonic tree are equivalent, except processing the image in 4-/8- (tree of shapes) or 6-connectivity (monotonic)
extrema-watershed tree: A. Vichik, R. Keshet, D. Malah: Self-dual morphology on tree semilattices and applications
constrained connectivity, alpha-trees, ultrametric waterseads: P. Soille, L. Najman: On morphological hierarchical representations for image processing and spatial data clustering


Can Oversampling help to mitigate the effect of additive white gaussian noise? Oversampling the signal at the receiver's ADC is a common technique to improve the ADC's resolution, dynamic range and to improve the signal to noise ratio (SNR) due to ADC quantizaiton noise, thermal noise,...
However, I am wondering if Oversampling in the ADC can help to reduce the effect of the additive white gaussian (AWGN) channel?
To my knowledge, oversampling in the ADC can help to improve the performance of ADC due to reduced ADC-induced noise, but I don't think it can reduce the effect of AWGN channel?
thanks

Oversampling will not do anything to change the existing noise in the waveform other than the quantization noise that is added. Sampling does not change the waveform but the quantization process of sampling with ADC’s will add additional noise to it (quantization noise and other additive noise collectively given as the “Equivalent Number of Bits” or ENOB). The OP is referring to the Gaussian noise that may be sampled, which itself will not be modified other than the additional noise added. This makes complete sense since our goal with sampling a waveform is to represent it with minimum distortion; there is nothing different between noise and a waveform of interest that would somehow be sampled differently by the ADC.
With regards to the quantization noise, oversampling will spread the additive quantization noise out over a wider frequency and thus reducing it density on a per Hz basis, since the quantization noise is well modelled as a white noise over the sampling bandwidth. The total quantization noise power doesn't change, but the noise density would be lower since that noise is spread over a wider bandwidth. When the ENOB is limited by spurs, then further oversampling will not reduce the noise further since a narrow band spur is confined in frequency rather than a spread noise. Further as sampling is increased the ENOB typically reduces due to non-linearities and aperture uncertainty in the ADC front end.
Given this, the typical goal in a system is to have sufficient gain in front of the ADC such that the amplified thermal noise floor is higher than the added noise from the ADC with sufficient margin such that the noise floor we observe digitally is the thermal noise floor of the analog waveform and not the added quantization noise of the ADC (for the latter case we would then be swamping out the signal with additional added noise which is typically not desirable when our interest is SNR). Oversampling can help meet this goal by effectively increasing the number of bits in the ADC (increasing its dynamic range).
For further details please refer to: What are advantages of having higher sampling rate of a signal?


May a causal system depend on future values of its own output For example, let's assume a system given as

y[n]-9y[n+2]=x[n]

I convert the difference equation to

y[n]=(y[n-2]-x[n-2])/9

and say that "the system is causal" and "not-memoryless" since it depends on past values of its output. But could it be possible to define a system to be "not-causal" and depend on any future value of its own output at the same time (no input term other than any present value of x[n] will be possible for this system and)?

By definition, the output of a causal system cannot depend on any future input or output values. However, I think there is a misunderstanding in your question. Suppose you are given the following difference equation:

$$a_0y[n+1]+a_1y[n]+a_2y[n-1]=b_0x[n]\tag{1}$$

The important point is that it is not possible to say if the system described by Eq. $(1)$ is causal or not, simply because Eq. $(1)$ does not uniquely describe one system, but in general there are three systems described by $(1)$. One of them is causal and its output can be computed from

$$a_0y[n+1]=b_0x[n]-a_1y[n]-a_2y[n-1]\tag{2}$$

Note that for any $n$ the output in $(2)$ only depends on past output (and input) values.

The second system described by Eq. $(1)$ is non-causal. The current output sample depends on past and on future output values:

$$a_1y[n]=b_0x[n]-a_0y[n+1]-a_2y[n-1]\tag{3}$$

And, finally, Eq. $(1)$ can also describe an anti-causal system, the output of which only depends on future output and input values:

$$a_2y[n-1]=b_0x[n]-a_0y[n+1]-a_1y[n]\tag{4}$$

The ambiguity of Eq. $(1)$ is also reflected in the $\mathcal{Z}$-transform domain, where $(1)$ corresponds to the following transfer function:

$$H(z)=\frac{b_0}{a_0z+a_1+a_2z^{-1}}\tag{5}$$

This transfer function, just like the difference equation $(1)$, does not uniquely define a system. Only if you specify the region of convergence (ROC) of $(5)$ will the system be uniquely characterized. Each of the three possible interpretations of Eq. $(1)$ discussed above correspond to a different ROC of $H(z)$.


Trying to find an intuitive understanding of OQAM/FBMC I am trying to understand, possibly in an intuitive way,  why would OQAM/FBMC when used with a Phydyas filter lead to a signal that does not need any cyclic prefix just like in any OFDM signal also created by an IFFT
If the cyclic prefix is absent how would ISI be avoided?

Since FBMC doesn't depend on the channel acting as circular convolution on the transmit signal like OFDM does, it makes no sense to use a cyclic prefix. (By the way, not all OFDM is cyclic-prefix OFDM, some OFDM-schemes simply use a silent guard interval, but these are the exception.)
To battle ISI without usage of a guard interval, one has to fall back to doing an equalizer per sub-channel in FBMC (something that OFDM reduces to multiplication with the inverse of the complex channel coefficient, but bought at the expense of rate loss due to necessity of the guard interval/cyclic prefix).
However, that equalizer is typically relatively benign in complexity, since the FBMC subchannels are designed to not see very many symbols in delay spread.
That's exactly what the PHYDYAS-style FBMC does: it uses the techniques used in modern single-carrier communications per carrier, and uses OQAM to combat inter-carrier-interference, whilst employing FFT methods for the multicarrier modulation.


Removing halo effect in image composition: Google Meet's solution Wanted to understand how google meet is removing halo effect when composing their image. They describe their method briefly in their blog as:

"The blur shader simulates a bokeh effect by adjusting the blur strength at each pixel proportionally to the segmentation mask values, similar to the circle-of-confusion (CoC) in optics. Pixels are weighted by their CoC radii, so that foreground pixels will not bleed into the background. We implemented separable filters for the weighted blur, instead of the popular Gaussian pyramid, as it removes halo artifacts surrounding the person. The blur is performed at a low resolution for efficiency, and blended with the input frame at the original resolution."

https://ai.googleblog.com/2020/10/background-features-in-google-meet.html

Am facing a similar issue where background color from original image is coming into the new composition as shown below creating a halo effect:

Am familar with guassian pyramids but not sure how it helping in this case and how seperable filters have an impact. Can anyone explain how to solve this problem?

Based on what I read in your reference I would say that
The key part for this solution is the model to compute the mask, if you have one that is a good starting point.
The second important feature is that it uses a position dependent filter, and the filter is parameterized with the length of the value of the mask. In their word avoiding the foreground to bleed in the foreground. They don't know how exactly, but I would say that from the mask a parameter $r_{i,j}$, that is large for the foreground and small for the background.
The possible realization of separable filter would be
$$ v_{i,j} = \frac{ \sum_k f(k\cdot r_{i,j}) x_{i-k,j}}{\sum_k f(k\cdot r_{i,j})} $$
$$ y_{i,j} = \frac{ \sum_k f(k\cdot r_{i,j}) v_{i,j-k}}{\sum_k f(k\cdot r_{i,j})} $$
And you are free to play with the normalization factors.
Using separable filter means that you use two 1D filters rather than small 2D filters for the image in multiple resolutions. This may be better in terms of memory. It can be implemented with two copies of the image, $x$ is written to $v$, and the second equation update the image $x_{i,j}$ in place.
Other possible advantage, since the filter is computed for every pixel separately. A separable filter of size $N \times M$, would require $N + M$ evaluations of the filter, while a 2D filter would require $M \times N$ evaluations, and in the realization above since the horizontal and vertical filters are the same it will require only $N$ evaluations.
Having the mask dependent filter is important to pass the foreground (the faces) with little distortion.


Does Reed-Muller codes exist when $r = 0$ and $m = 0$? Does Reed-Muller codes $\mathrm{RM}(r,m)$ exist for $r = 0$ and $m = 0$? Namely,  does $\mathrm{RM}(0,0)$ exist?  Some books mention that $m$ should be a natural number whereas some books mention that it can be a whole number including $0$,  for example Shu Lin and  Daniel Costello  Error Control Coding  book.

The codewords in a Reed-Muller code of degree $d$ and length $2^n$ are binary vectors of the form
$$\big(f(0,0,\ldots, 0, 0), f(0,0, \ldots, 0, 1),
f(0,0, \ldots, 1, 0), f(0,0, \ldots, 1, 1), \cdots, f(1,1,\ldots, 1,1)
\big)$$
where $f$ is the corresponding binary polynomial of degree at most
$d$ in $n$ variables.

Suppose first that $n > 0$. The dimension of the code is
$$k = \sum_{i=0}^d \binom{n}{i}\tag{1}$$
If $d=0$, then $k = \binom{n}{0}=1.$ There are only two "polynomials"
of degree $0$ in $n$ variables, namely the constants $0$ and $1$, and consequently
the codewords of length $2^n$ are $000\cdots 0$ and $111\cdots 1$.

If $n=0$ also, the codewords (if any) are of length $2^0 = 1$.
Now, you might want to claim from $(1)$ that $k = \binom{0}{0} = \frac{0!}{0!(0-0)!} = 1$ and so there are $2^k = 2$ codewords
(which, of course, are $0$ and $1$).
Alternatively, the $k$ in Equation $(1)$ is the total number of
subsets of cardinality $d$ or less of a set of $n$ elements.
There is exactly one subset (of cardinality $0$)
of the set of $0$ elements (namely,the empty set $\emptyset$).
and that subset is $\emptyset$ itself. Thus, the Reed-Muller code
of order $0$ and length $2^0$ is the $(1,1)$ code consisting
of two codewords $0$ and $1$. As a linear code, this is the
identity map.


Scatter Plot of OQPSK Below is my MATLAB code for OQPSK communication:

    data_in = [ 2 1 3 1 2 3 1 1 0 2 0 2 3 2 3 1 1 2 0 3 1 3  3 0 2 1 ];
    data_modualted=oqpskmod(data_in);
    filtmsg = filter(H,1,data_modualted); % Introduce channel distortion.
    eq1 = lineareq(nweigth,normlms(StepSize)); % Create an equalizer object. eq1.SigConst=oqpskmod(0:3);
    [symbolest,filtmsg_eq] = equalize(eq1,filtmsg,data_modualted(1:trainlen)); % Equalize.
    %%Scatter plot
    h = scatterplot(filtmsg,4,0,'g.');
    hold on ;
    scatterplot(symbolest,4,0,'kx',h);
    scatterplot(eq1.SigConst,1,0,'r*',h);
    legend('Filtered signal','Equalized signal','Ideal signalconstellation')

I am unable to understand which values to use to generate the phase
constellation diagram for my code?  How can i generate the correct
phase constellation diagram?

The first figure is constellation diagram of The filtered signal and the second is the constellation diagram of The Equalized signal

You'll need to plot your signal's imaginary and real (I and Q) part at the exact symbol time – in other words, you'll need timing recovery to get the constellation from a time signal.

Timing recovery is among the more complicated things you'll need to implement for a receiver, depending on the actual kind of transmission you're dealing with.

Assuming you weren't aware of that, you might want to  look into these terms (I think the order makes sense, didactically):

pulse shaping filter
eye diagram
Nyquist ISI criteria
timing recovery

square timing recovery
Gardner timing recovery

phase error correction

Costas loop


Difference between single sided and double sided amplitude spectrum I've recently gotten stuck with a small issue. I was asked to plot the single sided amplitude spectrum of the signal $x(t)=\cos(2\pi\cdot20000t)$.

I know that $X(f) = \frac{1}{2}\cdot\left[\delta(f-20000) + \delta(f+20000) \right]$. So there should be 2 impulses with amplitudes of $\frac{1}{2}$ at $\pm 20\,\text{kHz}$.

But I tried using Keysight's Advanced Design System to plot the single sided amplitude spectrum, and only got one impulse at $20\,\text{kHz}$ with an amplitude of $1$.

Why is this so?

For baseband real signals the frequency spectrum is conjugate symmetric; i.e., $$X(f) = X(-f)^*$$ This translates to an even magnitude spectrum; i.e., $$|X(f)| = |X(-f)|$$ and an odd phase spectrum.

Because of the symmetry it's also called as the double side band spectrum in commnnications terminology. For the real basband signal, a single side band is also sufficient to represent it. THere are different approaches to obtain the single side band as the uppersideband and the lower side band.

One method of obtaining the upper side band is to convert the signal into the analytic signal which is $x_+(t) = x(t) + j \hat{x}(t)$ where $\hat{x}(t)$ is the Hilbert transform of $x(t)$. The resulting spectrum is:
$$ X_+(f) =
\begin{cases}
2 X(f) ~~, &\text{ for} ~~ f > 0 \\
0 ~~, &\text{ for} ~~ f < 0 \\
\end{cases}
$$

Therefore for the given signal $x(t) = \cos(2\pi 20000 t)$ whose spectrum is $$X(f) = 0.5 \delta(f-20000) + 0.5 \delta(f+20000)$$ the upper side band will be given by $$X_+(f) = \delta(f-20000)$$

Note that you can adjust the amplitude of single side band if you wish to do so.


Having converted a grayscale image to binary using `im2bw()`, is the converse operation possible? I have converted a gray-scale image to binary using im2bw() now i wanna do the opposite. How can I do that in MATLAB?

As you've stated the problem, you can't. Converting a grey-scale image to black and white (binary) you're comparing pixel values (0-256 for an 8-bit image) to a threshold (that you define) and giving each pixel a value of 1 if the pixel is above the threshold, and 0 if it's below the threshold. The binary image doesn't contain information about by how much the pixel value was above or below the threshold.
If you want to get the pixel values from the original image which are above the threshold, you can multiply the binary image and the original image. All the values in the original image below the threshold will be zeroed out when they are multiplied by the zeros in the binary image. All the pixel values above the threshold in the original value will be multiplied by one and will retain their original values.
Unfortunately, if you're given an arbitrary binary image with no information about the threshold used, or the original image, you cannot recover the original grey-scale values.
Hope this helps.


Viterbi decoder on a bit stream where the preamble is also coded I hope this is not off-topic for this site.

In my problem, Viterbi decoder is used for decoding convolutional-coded (rate = 1/2, constraint length = 7) messages.

A message is 250 bits where the first 8 bits are 1/3 of a preamble. There are several types of messages. A message is received as 500 bits because of the coding, which means, the preamble is also coded.

Since there is no place marker other than the coded preamble, I am in trouble with finding the first bit of the 500 bits to start decoding.

I would really appreciate any suggestions and comments.

Is the 500-bit encoded message part of a continuous stream being decoded by the Viterbi decoder, or is it a 500-bit packet that is decoded individually?   If part of a stream, then decoding continuously and searching for the preamble indicates the beginning of the message.   If packetized, then the decoder memory needs to be managed so that the first decoded bit location is known.   This will be dependent on the decoder implementation.