Filter Design near Nyquist frequency (complex curve fitting) I'm fitting IIR filters to a complex transfer function with Matlab's invfreqz.m. My fitting results are very good, but when I come very close to the Nyquist frequency the phase response runs to 0 and I also get deviations in the magnitude response. I did some research but couldn't find anything regarding that phenomenon. Can anybody point me in the right direction, for which keywords I have to search in order to find something on that topic and how I can try to optimize the fitting near the Nyquist frequency?

You can't get a good fit very near the Nyquist frequency (or near DC) because information is lost due to the imaginary components around that frequency being of (almost) opposite phase and thus (nearly) cancelling out.  Thus, there will be ambiguity errors due to trying to fit any curve to lost information.

If possible, you can try to gather more information near (but not at) the Nyquist frequency by sampling for a longer time span or at a higher sample rate.  This will move the measurement point slightly farther from where there is complete cancellation of the imaginary component.


attenuation of the central difference differentiator In the following snippet, I am differentiating a sine wave using the central difference equation,
first through misc.derivative and then through convolution with the kernel [1/2, 0, -1/2]
import numpy as np
import scipy.misc as misc
import scipy.signal as sig
ncycles = 5.0
nsamples = 500.0
x = np.linspace(0,ncycles*2*np.pi,nsamples)
dydxmisc = misc.derivative(np.sin,x,dx=0.1)
dydxconv = sig.convolve(np.sin(x), [0.5,0.0,-0.5], mode='same')

As shown in the plots below, misc.derivative produces the expected output when dx is small enough, but when using convolution the derivative is attenuated.

The attenuation changes based on number of cycles and number of samples and is a factor of pi but I cannot figure out what the exact relationship is. Can anyone explain why the attenuation happens and what the relationship is? In other words, how do I have to scale the kernel in order to get a derivative of the correct dimension?

Your filter forgets that in the discrete differences equation, you divide the difference in value by the difference in $x$; so, you'd need to multiply the result with ncycles*2*np.pi/nsamples to get the same.
Other than that, you're not compensating for the group delay of your FIR filter, either, so you want to shift the result by half the filter length, too.


How does nearest neighbour, bilinear and cubic interpolation work in images? More math is appreciated for each of the methods and references are appreciated. I have tried understanding from Wiki and matlab link but don't understand how the translation matrix is being used there.

What I seek is to be able to teach this to someone else. Thanks!

Here is an illustration of bicubic interpolation, which typically needs 4x4 input pixels to calculate one output pixel:

Figure 1. Bicubic interpolation from 4x4 input pixel values to calculate one output pixel value. a) The output pixel must be located in the greyed area with respect to the input pixels (black dots). b) Four intermediate values (red dots) are independently calculated by interpolating horizontally from each row of four input pixels. c) Using the same interpolation method, the output value (blue dot) is interpolated vertically from the column of four intermediate values.

Each interpolation amounts to a weighted sum of four values, with the weights given by cubic polynomials of horizontal or vertical position within the grey square.

Bilinear interpolation works similarly but only has as input the 2x2 pixels at the corners of the grey square, and the polynomials are linear.

Nearest neighbor interpolation has the grey square centered at a pixel, and simply that pixel value is output.


Proof of Paley-Wiener criterion for causality The Paley-Wiener criterion for causality is that $\displaystyle\int_{\mathbb{R}}\frac{A(\omega)}{1 + \omega^2}\mathrm{d}\omega$ exists and is finite, where $A(\omega) = \left|\mathcal{F}[f]\right|$ is the modulus of the Fourier transform of a function $f: \mathbb{R}\to\mathbb{R}$

How is this derived from the usual statement of the Paley-Wiener theorem?

I have also seen that this proves that a causal filter can't provide total rejection for a band of frequencies, but this appears to me to follow directly from the usual statement that the Fourier transform of a square-integrable distribution vanishing for negative reals is holomorphic and the identity theorem for holomorphic functions.

This is nice question which I have recently encountered while studying ideal filters. Suppose transfer function $H(p)$ which is complex rational and holomorphic (or meromorphic?) in $p = \sigma + j\omega$ and satisfy Cauchy-Riemann conditions.

We can rewrite transfer function as $H(p) = e^{-\Gamma(p)}=10^{-\frac{1}{20}\Gamma(p)}$ with respect to nepers and decibels. Since exponetital function is holomorphic so the $\Gamma$ is. If we restrict ourselves to imaginary axis $p = j\omega$ we have

$\Gamma(j\omega) = A(j\omega) + jB(j\omega)$

where $A(j\omega)$ is an attentuation charateristics and $B(j\omega)$ is a phase characteristics.
In general we have identity $\ln z = \ln |z| + j\arg z$ for any $z \in \mathbb{C}$. So we can write

$A(j\omega) = \ln|H(j\omega)|$ and $B(j\omega)=\arg H(j\omega)=\arctan\frac{\Im\lbrace H(j\omega)\rbrace}{\Re\lbrace H(j\omega)\rbrace} $.

Then consider pair of Hilbert transforms given by integrals (Cauchy principal value are taken in each of the integrals)

\begin{equation}
A(\omega) = \frac{1}{\pi}\int_{-\infty}^{\infty}\frac{B(j\xi)}{\omega - \xi}\,\mathrm{d}\xi \quad\quad B(\omega) = -\frac{1}{\pi}\int_{-\infty}^{\infty}\frac{A(j\xi)}{\omega - \xi}\,\mathrm{d}\xi
\end{equation}

Evaluating these integrals for module characteristic of ideal low pass filter - say let $A(\omega)$ is constant at $a$ from $-\omega_1$ to $\omega_1$ (module characteristics are symetric) gives
\begin{equation}
B(\omega) = - \frac{a}{\pi}\ln\frac{|\omega + \omega_1|}{|\omega - \omega_1|}
\end{equation}
Thus $B(\omega)$ has points in infinity at $\pm\omega_1$.

See picture bellow:

Similarly for linear phase characteristics $B(\omega) = a\omega$.

This is saying to us that we can't reach filter transfer function $H(p)$ with both ideally rectangular module characteristics and ideally linear phase characteristics at the same time.

We see that $B(\omega)$ with points at $\pm \infty$ can't realise any phase characteristics $\arg H(j\omega)$ of a given feasible filter. That's the reason why
\begin{equation}
-\frac{1}{\pi}\int_{-\infty}^{\infty}\frac{A(j\xi)}{\omega - \xi}\,\mathrm{d}\xi < \infty
\end{equation}
must converge and so the
\begin{equation}
-\frac{1}{\pi}\int_{-\infty}^{\infty}\frac{|A(j\xi)|}{\omega^2 + \xi^2}\,\mathrm{d}\xi < \infty
\end{equation}
must.


Converting ADC output to fixed-point format I want to convert the two's complement, 16-bit output from an ADC into a generic fixed-point format, $\rm Qm.n$, where $\rm m$ is the number of integer bits and $\rm n$ is the number of fraction bits.

Any hints/procedure to do this job would be appreciated.

You don't have anything to do, because $\rm Qm.n$ is just an "interpretation" of an integer. The only thing that really matters is that you need to correctly apply a shift operation when multiplying numbers; or when adding/subtracting numbers with different values of $\rm (m, n)$ - You can think of the $\rm Qm.n$ representation as a "tag" you attach to number, a bit like physical units.

In your case, the output of your ADC is an integer between $-32768$ and $32767$. This number can be interpreted as a $\rm Q15$ fixed-point number between $-1.0$ and $0.99997$. Or as a $\rm Q2.14$ fixed-point number between $-4.0$ and $3.99987$. Or as a $\rm Q8.8$ number between $-256.0$ and $255.992$. You don't have to do any arithmetic operation.


FFT of long signal by segments/chunks with discontinuities I am processing an EEG brain signal, which has up to 64 data channels sampled at 500 Hz. One of the analyses consists of extracting the ratio of alpha/delta power, where alpha represents the waveforms which frequencies are between 8 and 13 Hz, and delta the waveforms which frequencies are between 1 and 4 Hz.
The signal is first filtered with a 1-15 Hz BP filter (scipy Butterworth 4th order, output 'sos'). Additionally, a common average projector and EOG SSP artifact correction projector are applied. Those are not relevant to this question.
My program compute the alpha/delta ratio at 2 different moments:

A short window (500 ms or 1 second, i.e. 250 or 500 samples)
A long window (4 to 10 seconds)

However, EEG signal is very susceptible to noise, i.e. if the participant moves his eyes, jaws, head, .. the brain signal will be completely masked by artifacts. Thus, a rejection criteria is also applied. It's a simple peak-to-peak rejection, i.e. if the max - min is larger than the criteria on at least one of the channels, this window (epoch) is rejected.
For the short window, I can directly implement this as follows:
if any(np.ptp(data, axis=0) > reject['eeg']):
    continue # skip

With data a channel x samples array (e.g. 64x250) and reject['eeg'] the threshold value.
However, for the large window, I can not disregard the complete window. My idea was to cut down the e.g. 5 seconds signal into 1 second chunks, and apply the rejection criteria on the chunk.
Something like:
data = raw.get_data(picks='eeg') # np.ndarray 64x2500 for 5 seconds
data = data.reshape(64, 5, 500) # I hope this is the correct reshape, not tested.
chunk2remove = np.any(np.ptp(data, axis=2) > reject['eeg'], axis=0)
data = data[:, np.where(chunk2remove == False)[0], :]

However, now I do not know if I can safely apply FFT on the signal. I do not know what the impact of this chunk rejection is, nor what the border effect might be. What is the correct way to apply an FFT to this chunked signal with discontinuities between the chunks?
Currently, the method I use to extract the alpha/delta power ratio is:
fs = 500.
window = 5 # window length in seconds
fft_freq = np.fft.rfftfreq(int(fs * window), 1.0/fs)
alpha_band = np.where(np.logical_and(fft_freq>=8, fft_freq<=13))[0]
delta_band = np.where(np.logical_and(fft_freq>=1, fft_freq<=4))[0]

fftval = np.abs(np.fft.rfft(data, axis=1) / fs)
alpha = np.average(np.multiply(np.abs(fftval[:, alpha_band]).T, weights))
delta = np.average(np.multiply(np.abs(fftval[:, delta_band]).T, weights))

The array weights is simply a [0, 1] weight applied on each channel. As you can see, this is not exactly the power, but as I'm doing a ratio, I did not feel that the squaring was required. Feel free to correct me if I'm wrong.
If you have any comment, additional ideas for improvement, I'm also interested.

I can think of two things to consider: discontinuities and phase.
The first is that whenever you take a finite chunk from some conceptually infinite signal, you get a step-discontinuity (ie. click) in the envelope of the signal, which generally results in some broadband energy. This can be reduced by multiplying the chunk with a window function, at the cost of slight blurring of the desired spectrum). I don't see any windowing in your code, so you might want to at least think about this. Standard windows might not be ideal though, see below.
The other thing to consider is that if you take a larger chunk and remove a part of it form the middle, then there will again be a click. If the two chunks are moved next to each other, then the phases from the two chunks probably won't align very well and the measured magnitudes will be lower in a frequency dependent way (probably unwanted, since this is generally "unpredictable").
That said, the phase-problems at least can be avoided if rather than removing a chunk, you simply set it to all zero instead so the time-alignment is maintained (ie. temporarily mute the signal; I apologize if you are already doing this, my Python is rusty). The muted part obviously won't contribute, but otherwise the only remaining impact I can think of should be the extra envelope transients.
To reduce the impact of the envelope transients (assuming the broadband energy causes problems), you might want to do some fade-in/fade-out, for example by using the cosine-tapers of a Tukey-window wherever it starts or ends (ie. around muted parts + beginning and end of the whole thing), though I don't understand your specific application well enough to really give more specific advice than this.
Finally I want to note about the squaring: while $a^2/b^2 = (a/b)^2$ the average of squared magnitudes (or the root of average of squared magnitudes) is not the same as the average of absolute magnitudes; the relative weighting of values is different. If you want to compute the average power, then you should square before averaging.


How i can create ghost on wavelet or seismic trace and then Deghosting that in MATLAB? I want to model a seismic trace or wavelet and then create ghost on them, at the end Deghosting in MATLAB.

Can anyone help me please?

YOu can use SeisLab, which apparently has an option for that:

  0,6.1,'s_ghost_response','Apply effect of source/receiver ghost on seismic dataset';


How do I initialize the state of a digital filter in Direct Form II? Suppose I have a digital filter implemented in Direct Form II. How do I initialize the state of the filter as if the input $x[n]$ had a fixed value $x_0$ for all $n<0$?

The difference equations for this filter are:

y[n] = b0 w[n] + b1 w[n-1] + b2 w[n-2]
w[n] =    x[n] - a1 w[n-1] - a2 w[n-2]

To achieve steady-state, w[n] == w[n-1] == w[n-2]. Call this value w. Solving the second difference equation, we find w = x/(1 + a1 + a2), where x is the steady-state input value.


Calculate the time delay introduced by group delay for IIR-Filters First of all, let me show what I mean by a simple example of an FIR-Filter:

I designed an FIR-Filter with a constant group delay of 323 samples and my sampling rate is 500 Hz. Now I can calculate the time delay by this filter with 323/(500 Hz) = 0.646s.

An IIR notch filter has a frequency dependent group delay and I am interested in the time delay around 50 Hz.

We know the group delay is around 2.5 samples at 50 Hz and we sample with 500 Hz. Can I do the same calculation as before? This means the time delay for 50 Hz is 2.5 samples/500 Hz = 5ms.

I am not sure but I think the correct calculation is with fs/2, that means it's 2.5/250 = 10ms.

If you want to compute the delay experienced by a sine wave at a given frequency $\omega_0$, then it's not the group delay but the phase delay that you need to consider. In the case of linear phase FIR filters the phase delay and the group delay are identical, but in the general case they aren't.

The phase delay is defined by

$$\tau_p(\omega)=-\frac{\phi(\omega)}{\omega}\tag{1}$$

where $\phi(\omega)$ is the frequency dependent (unwrapped) phase of the system's frequency response. Clearly, if $\phi(\omega)=c\cdot \omega$ with some constant $c$, then the phase delay $(1)$ equals the negative derivative of the phase, i.e., group delay and phase delay are equal.

Of course, the (de-)normalization by the sampling frequency is always done in the same way. If $\omega$ is normalized, then you get the phase delay in samples, otherwise in seconds.

EDIT:

The response of an LTI system with frequency response $H(\omega)=M(\omega)e^{j\phi(\omega)}$ to a sinusoidal input signal $x(t)=A\sin(\omega_0t+\theta)$ is given by

$$\begin{align}y(t)&=AM(\omega_0)\sin(\omega_0t+\theta+\phi(\omega_0))\\&=AM(\omega_0)\sin\left[\omega_0\left(t+\frac{\phi(\omega_0)}{\omega_0}\right)+\theta\right]\end{align}\tag{2}$$

So the delay experienced by $x(t)$ is given by $\tau_p(\omega_0)$ as defined in $(1)$.


Breaking a convolution into smaller pieces For a project I need to do convolution and i use gpu for calculations. Sometimes I have to deal with kernel sizes of 50x50 and this size of kernel is sufficiently large that it chokes the gpu. (not enough memory svailable by the gpu) I need to find a way to break the kernel into smaller sizes (8x8 or similar) and do the convolution this way (i.e. piece by piece and later stitch together) so that I can realize gpu enhancements. What is the best way to break these kernels into smaller pieces? (I have no prior knowledge of kernel size so the solution must be something I can deal with in run time. )

Yes, you can split them up.  Convolution is a linear process, which means that superposition holds.  Thus, you can break up any convolution kernel $k$ into multiple parts ($k_1, k_2, ... k_N$) such that $k = \sum k_i$.

For example, if you had a convolution kernel that looked like [1 2 3 4] (granted this is a silly kernel, but does fine for purposes of illustration), you could break it up into the following kernels, [1 2 0 0] and [0 0 3 4].  Now, you could simply do the convolutions like normal and then add the results together.  The sum would be equal to the convolution product of the original kernel.

Doing that would be inefficient, though, both in terms of computations and memory.  Instead, you could simply drop the 0's at the ends of the "sub-kernels", and do the convolutions using [1 2] and [3 4].  The one tricky part to this is that you have to offset the results before adding them together to get the correct answer.


Deciding if the received packet is valid or not valid I have sampled RF wireless packets as IQ samples file, those RF packets has patterns like preamable|syncword|data|crc , I'm doing cross correlation with my constant syncword (after I did the demodulator in order to convert the IQ samples ) and I see in matlab the peak of correlation where there's sync word in my packets, Im now trying to decode those RF packets which I want to decide if the packets that I received are valid or not valid.
what should I do in order to decide if the packet is valid or not?? maybe it's related to CRC but how can I decide if the packet is valid or not? what I have done and arrived to this stage : I did cross correlation with my syncword and I see the peak of corelation of syncwords !
Anyone could help me by the flowchart what should I do in order to decide if my received packet are valid or not valid? thanks
CRC/SYNCWORDS in the transmitter are constants ..

CRC is a "definite" indicator of the channel impacting the bits and error in received symbol. So deploying a CRC check would definitely help. Length of CRC would depend on the overall system and data packets.


Extremes of integral of the autocorrelation theorem I am trying to apply the correlation theorem (here attached) to this definition of correlation of a finite continuous signal $g(t)$:
$\rho(\tau)=\frac{1}{T} \int_{t_{0}}^{ t_{0}+T}g(t)^{.}g(t+\tau)dt  \tag{1}$
where $t$ is the time, T is the time duration of the signal, $t_{0}$ is the starting time, $\tau$ is the time shift of the correlation.
As you can see, in the attached demonstration the correlation is defined from $[-\infty,+\infty]$. Obviously, the signals we work with are always finite. So I defined the interval as $[t_{0},t_{0}+T]$, as I also see in a book of Turbulence.
If I do the fast Fourier transform of both the member of the Eq (1) I obtain
$E(f)= \frac{1}{T} \int_{\color{red}{t_{0}}}^{\color{red}{t_{0}+T}}\int_{\color{red}{t_{0}}}^{ \color{red}{t_{0}+T}}g(t)^{.}g(t+\tau) e^{-i2\pi f \tau}d\tau dt \tag{2}$
where $E(f)$ is the Fourier Transform of $\rho(\tau)$. Usually the Fourier transforms, that I have seen in books, are defined in the interval $[-\infty,+\infty]$, but, as I said before, the signal we work with are defined in a time window such as $[t_{0},t_{0}+T]$, so I tried to apply them to finite signal
Introducing the variable
$y=t+\tau$
The Eq. (2) can be rewritten:
$E(f)= \frac{1}{T} \int_{\color{red}{t_{0}}}^{\color{red}{ t_{0}+T}}\int_{\color{red}{t_{0}}}^{ \color{red}{t_{0}+T}}g(t)^{.}g(y) e^{-i2\pi f (y-t)}dy dt \tag{3}$
Furthermore the Eq. (3) can be rewritten as:
$E(f)= \frac{1}{T} \int_{\color{red}{t_{0}}}^{\color{red}{ t_{0}+T}}g(t) e^{+i2\pi f t}dt\int_{\color{red}{t_{0}}}^{ \color{red}{t_{0}+T}}g(y) e^{-i2\pi f y}dy \tag{4}$
So E(f) can be rewritten as:
$E(f)= \frac{1}{T} G^{*}(f)^{.}G(f)\tag{5}$
I am not sure about the (red) extremes of the integral, I think they are wrong can you help me? Is it possible to apply the correlation and the Fourier transform to finite signals with interval of existence $[t_{0},t_{0}+T]$?

Maybe your signal is finite, but the function you are transforming is infinite. So you can do your demonstration using the interval $[0,T]$ or $[-\infty,+\infty]$, the result should be the same, because if your signal is finite the FT threats it as an infinite signal with period T.


How to draw a finite-duration sequence On page $359$ of the 2nd edition of the book "Signals & Systems", by Oppenheim and Willsky, the authors wrote the following words:

"Consider a general sequence $x[n]$ that is of finite duration. That is, for some integers $N1$ and $N2$, $x[n] = 0$ outside the range $-N1\leq n\leq N2$. A signal of this type is illustrated in $\text{Figure}\; 5.1\:(\text{a})$"

My three questions are:

Regarding $x[n]$ in the figure, do the leading and trailing zero-valued samples exist?
If the leading and trailing zero-valued samples exist, then what is the duration of (how many samples are in) the finite-duration $x[n]$ sequence?
If the leading and trailing zero-valued samples do not exist, then why are they shown in $\text{Figure}\; 5.1\:(\text{a})?$

The plot as shown makes sense as nothing has limited the time axis to be going from $-\infty$ to +$\infty$. So we are showing the finite duration sequence as it would appear on such a time axis.

We can limit the time axis to a finite number of samples as used in the Discrete Fourier Transform (DFT) or the an unlimited time axis as used in the Discrete Time Fourier Transform (DTFT). Neither approach would change the number of samples in a finite duration signal.

So the plot and text is implying the time axis extends from  $-\infty$ to +$\infty$. The number of samples in the finite duration sequence would be the non-zero samples.

My confusion when first answering this question was the case of the DFT and how adding zeros is significant to the result, but really this is just approximating the DTFT and not changing the duration of the sequence itself.


How to implement simple noise filter based on the frequency domain information? Brief introduction:
I am working on my university project. My task is to develop application that captures sound from mobile phone and removes/reduces noise (this simpilied version of task).

I never faced with DSP before, I spend a lot of time to figure out how to get audio signal from phone (no direct API is porvided), so I have only few weeks left to study DSP.
There are a lot of information on that field, my basic source of information was: http://www.dspguide.com/.

My work, to be accepted should include some work(algorithms) in the field of DSP (so I should avoid ready libraries).

What I need
simple algorithm(easy to understand) for noise reduction, I am very limited in time (I got only one day left, after that I should present my work)

refernces to resources for understanding FFT & inverse FFT (although I implemented that algorithms in code, I simply adopted code from the book)

what I have done
FFT, Inverse FFT
I tried to remove frequencies by applying FFT to signal in time domain and zeroing needed ranges in result frequency domain (ranges < 300 && ranges >3700 - as human voice is in range [300;3700] - so doing that I tried to remove sounds not related to human voice, but I got very bad results (https://stackoverflow.com/questions/24101814/why-ideal-band-pass-filter-not-working-as-expected)

Could you please suggest me how I can based on the results of FFT reduce noise (I don't need super efficient method, I need something that I can present in my work and explain how I did it)

Please help!

Fast convolution filtering with an FFT/IFFT requires zero padding and using overlap add or overlap save methods to remove circular convolution artifacts.

You should also use a frequency domain response that has a shorter impulse response than a rectangle (which zeroing bins is).


How would you generate this type of wave form? I'm trying to generate a wave form that looks like this:

It looks to me like a square wave that has a large transitional section between going from "low" to "high" and back.

Any idea how to generate such a wave mathematically (and/or via some Python library)?

Note also that there is a higher-frequency carrier wave.

Try a suitably parameterized function of form:

$$f(t) = \operatorname{clip}\big(a_0\sin(f_0 t) + a_1\sin(f_1 t)\big),$$

where $t$ is time, $f_0$ and $f_1$ are the frequencies and $a_0$ and $a_1$ the amplitudes of two sine waves, and:

$$\operatorname{clip}(x) = \left\{\begin{array}{ll}-1 &\text{if } x < -1\\x &\text{if }-1 \le x \le 1 \\1 &\text{if } x > 1.\end{array}\right.$$

You don't show the scale of the data so the bounds $-1$ and $1$ might be of wrong absolute value. I assumed that the data is zero-centered so there is no zero frequency (bias) term and the bounds look symmetrical in your graph.


Application of MAP criteria This question has been asked in GATE 2019. I am unable to understand how MAP criteria is "applied" here. I've read about MAP but I'm unable to get how they're solving it by graphical method.
Please explain what they are doing here, and also what do they mean by "optimum threshold exists here" ?

Thank you.

The MAP criterion minimizes the average probability of error, which in your case reads
\begin{equation}
\begin{aligned}
P_e &= P \left( -1 \right) P \left( Y>\theta \mid -1 \right) +
P \left( +1 \right) P \left( Y<\theta \mid +1 \right) \\
&= 0.2 \int_{\theta}^{+1} f_Y \left( y \mid -1 \right) \mathit{dy} +
0.8 \int_{-1}^{\theta} f_Y \left( y \mid +1 \right) \mathit{dy} \,.
\end{aligned}
\end{equation}
Therefore, the optimum threshold, call it $\hat{\theta}$, is optimum in terms of minimizing $P_e$.
Exploiting the previous equation and the definitions of the random variables $X$ and $N$, you get
\begin{equation}
\begin{split}
P_e &= 0.2 \times 0.25 \times \left(1-\theta\right) +
0.8 \times 0.25 \times \left(1+\theta\right) \\
&= \frac{1}{20} \times \left(1-\theta\right) +
\frac{4}{20} \times\left(1+\theta\right) .
\end{split}
\end{equation}
Minimizing $P_e$ w.r.t. $\theta$, you obtain $\hat{\theta} = -1$.
An easier way to solve your problem would be: look at the plot, and select $\hat{\theta}$ such that it minimizes the sum between the right tail of the left function (i.e., the grey area) and the left tail of the right function. In this way you easily get the solution, but the analytical way is, in my opinion, preferable to understand the concepts.


Stereo sound separation into three sources Whenever we listen stereo songs, one part of the audio like voice hear at the center of the head, and other parts are nearby left ear and right ear. If I represent center part by C, and other parts by L' and R', We may model our original left and right channel as,

L' + C = L    and,

R' + C = R.

As, common part of L and R channel comes in C vector, L' and R' are uncorrelated. Now we can say, (alpha + beta = pi/2)

Is this modelling correct as per the listening scenario for stereo sound?

Partially correct. Human auditory spatial perception is rather complicated. A good book on the topic would be Spatial Hearing.  How and where we localize things depends on a great many factors. For left/right localization the main effects are interaural time difference and interaural level difference. Time differences is the strongest queue. Whatever comes earlier on the left ear will be solidly localized on the left unless even if the sound is louder on the right (within reason).
If left & right signals are perfectly time aligned than your model is simplified model.


Signal Acquisition in Compressed Sensing I'm trying to wrap my head around compressed sensing, so I've been reading this intro to the topic.

I'm completely keeping up when they discuss exploiting the sparsity of a signal in some domain to compress it.  I also think I get why you'd want to choose a sensing matrix that is incoherent with the basis in which your signal is sparse.

The point where my understanding breaks down is the "Undersampling and Sparse Signal Recovery" section.

Concretely, let's say I have a signal f of length n.  I have a sensing matrix which is incoherent with the basis in which f is sparse that is m x n, and m < n.

If I understand correctly, I'm supposed to measure the signal by taking its inner product with each row of my sensing matrix, so I'll end up with samples of length m.

As I understand it, compressed sensing is useful in cases when sampling at the Nyquist/Shannon frequency is prohibitively expensive or slow.  I must be missing something really fundamental, but didn't I have to do the work to acquire all samples of f so that it could be correlated with my sensing matrix?  As I understand it, I've compressed the data now, which is great, but I didn't save any energy in the acquisition process.

Just to restate everything for clarity, the Compressed Sensing problem is defined as the following: given a signal $x$ of length $N$, we measure the projection of $x$ by some projection operator, $\Phi$ of size $M \times N$,

$$ y = \Phi x,$$

where $y$ is are our $M$ measurements of the original signal. We could also state this in terms of the inner products of $x$ with the rows of $\Phi$,

$$ y_i = \left< \phi_i, x\right>.$$

In the CS sampling procedure we have an inherent dimensionality reduction since we are taking a projection of $x$. The degree of this dimensionality reduction is usually referred to in terms of the ratio $\frac{M}{N}$ (aka the subsampling ratio or "Subrate"). By obtaining $y$ during acquisition we are simultaneously performing acquisition and dimensionality-reduction rather than full-resolution sampling followed by compression (DCT, DWT, etc.). The measurements $y$ are read off of the sensor, quantized (still an open problem, but decent results can be obtained through simple scalar quantization. See Laska's dissertation for more novel approaches), entropy coded (pick your flavor), and then transmitted or stored.

You touch on a good key point when it comes to CS signal acquisition, namely, if the calculation of this projection is costly then the advantage of CS for acquisition systems would appear small. However, depending on the context and type of signal, such projections can be accomplished in the analog domain (requiring no computation). The single-pixel camera (SPC) is a great example of this.

Researchers have also been able to greatly reduce sampling times for MRI using CS techniques. Specifically, conventional MRI techniques sample along different radial lines within the frequency domain. Each radial line is a measurement of the MRI device requiring some sampling time. Traditionally, for higher resolution MRI, more radial lines must be acquired, incurring longer acquisition time (problematic for the MRI of small children). However, this process can also be characterized as linear projections within the frequency domain. Because of this, CS techniques can be used to recover MRI images from many fewer measurements, allowing for high resolution MRI with much shorter MRI appointments. This was, in fact, the context in which CS was first applied. For more info, this paper by Lustig et al might be a good starting point.

In summation, the usefulness of CS to a particular context really depends on how you implement the projection. Getting it right can require some out-of-the-box thinking to come up with a novel sensing strategy. Thankfully, frameworks such as the SPC are generalizeable to a wide range of different signal contexts, so we don't have to re-invent the wheel every time.


Confusion regarding notation/symbol? I was reading circular convolution topic from Digital Signal Processing Using Matlab  3rd Ed,by Proakis. I came across a strange term/symbol; I have highlighted it in attached snapshot
Update:Pg150; same chapter,same book(just before example 5.4) contains this notation and mentions that it is a rectangular window of length N
I have also attached snap of eq 5.24,in 2nd snap, where this RN(n) symbol appears also:

This symbol denotes a rectangular pulse of length $N$:
$$\mathcal{R}_N(n)=\begin{cases}1,&0\le n\le N-1\\0,&\textrm{otherwise}\end{cases}$$
I'm not sure where it is defined for the first time, but this definition is clear from the equation above Eq. $(5.24)$ on page $130$ of the 3rd edition of Digital Signal Processing Using Matlab by V.K. Ingle and J.G. Proakis.


Why there is Difference between shapes of ROC of z domain and s domain? ROC(region of convergence) of Z domain is shown by a circular region while ROC in S domain is shown by a rectangular(approximately looking like rectangle) region

What is the reason of this difference in shapes of ROC regions?

Because, the region of convergence in the Laplace transform
$$
X(s) = \int_{-\infty}^{\infty} x(t) e^{-st} dt
$$
is related to the weighting provided by the real part of the complex $s = \sigma + j \omega$; as this will yield the weight $|e^{-st}| =  e^{-\sigma t}$  applied on the input signal $x(t)$, and is a function of $\sigma$ alone and is a rectangular (planar) region on the s-plane.

But the region of convergence in the Z-transform
$$
X(z) = \sum_{n=-\infty}^{\infty} x[n] z^{-n}
$$
is related to the weighting provided by the magnitude of the complex $z = \sigma + j \omega$ as given by $|z|^{-n} = |\sigma + j\omega|^{-n} = |z|^{-n} |e^{-j n \angle{z}}| = |z|^{-n}$, which is a circular region on th z-plane, due to the magnitude of $z$ being involved.


Carrier sense multiple access, MAC layer A medium access control MAC mechanism in wireless local area network and Ethernet is called Carrier Sense Multiple Access.

There are two different types of this CSMA

1- CSMA/Collision Detection

2- CSMA/ Collision Avoidance

Does anyone know what the advantages and disadvantges of each are?

In collistion detection, a node may start transmitting at any time, and backs off if it detects that it caused a collision. This forces the node that was already transmitting to back off, too. Its advantage is simplicity; its disadvantage is that it introduces big inefficiences: a packet makes it through the network only if no other node decides to transmit at the same time.

In collision avoidance, a node first looks at the medium. If it detects a signal, then it backs off and waits before attempting to transmit again. This increases the transmitter's complexity a bit, but the network efficiency is increased. A collision only happens when two nodes decide to transmit at (almost) the same time.

You can find the details and how to calculate the specifics in any textbook on networking.


Resolving two nanosecond pulses with filter Problem:

I have to design a filter that should be able to resolve two ~20ns pulses as well as get rid of some ringing. I was able to find a filter kernel that if convolved with the noise signal is able to do the required filtering. Below are the graphs of the pulses before and after applying the convolution.

The graph below is the filter kernel that I used to do the filter in both the time and frequency domain.

Attempted solution:

According to my understanding -- which is very superficial -- the discrete Fourier transform of the filter kernel is the frequency response of the filter. Based on that I went on and tried to design a filter with a similar frequency response which according to the graph is a low-pass filter. However, when I design and apply a low-pass iir filter to the signal it is not able to resolve the pulses. What's wrong with my understanding? How could I go about designing a filter that has a similar filter response as the one I require? Is it possible to do it with an analog rather than a digital filter?

PS: I was not able to upload more pictures but if possible I will upload them to the comments as well as the code and data use to produce the graphs.

The problem is that your low-pass IIR filter matches the original filter's magnitude response but will definitely not match its phase response.

The IIR filter will have a kernel that is larger closer to $t=0$ than your plotted filter kernel.  This means that, even if the magnitude frequency response matches, the phase frequency response will not be similar.

So you have the filter's impulse response, which is sort of like a long FIR filter.

Unfortunately, designing an analog filter that matches this response is a little tricky. Most techniques for designing such filters assume you don't care about the phase response.

One way might be to take your FIR current filter and perform model order reduction on it. See, for example, matlab's balred method that does balanced truncation (I think), though you may need to do some discrete-to-continuous or vice-versa transformation to get things to work.


Easiest Pattern to Recognize with Machine Vision I need to detect the exact position of some kitchen tools using machine vision. I've the possibility to label the tools with a pattern, to make the detection (more specifically the image segmentation) easier.
What is the easiest pattern to detect using machine vision?
The pattern should be:

Easy detect in presence of partial occlusions.
Easy to detect with varying and degraded light conditions.
Easy to detect with widely used computer vision algorithms, in particular lightweight algorithms that can run on edge devices.
Hard to mistake for another object (in a kitchen).
Possibly carrying some additional information.

A QR code fits this description pretty nicely. My guess is that the easiest pattern will contain a sharp alternation of two colors, likely black and white.

In case you can shoot a video of the static scene than a blinking light would be the easiest as you could easily detect it by subtracting the n - 1 frame from the n frame until you see something with high values.
If you take a still shot you can use 2 main ideas:

If the colors of the scene are from a given plate, find a color very different in Hue and make sure it is intense (Even active lighting which will saturated).
Use a pattern which is easily detectable like stripes of B/W or pure green and magenta. You may encode data in the stripes.

Resources
Some resources pasted here in comments:

Fiducials / Fiducial Markers.
AprilTag 3 - A visual fiducial system popular in robotics research.
What the Name of This Kind of Computer Vision? (Related question).


Polyphase Decomposition in Matlab I created the lowpass filter

h = firls(100, [0 .2 .25 1], [ 1 1 0 0],[.001 60]);

My goal is to get four filters  $H_{k}(z) = H(W_{4}^{k}z)  k =0,1,2,3$

From my knowledge to get the polyphase components you would do the following

$H_{0}(z) = \sum h(4n)z^{-n}$
$H_{1}(z) = \sum h(4n+1)z^{-n}$
$H_{2}(z) = \sum h(4n+2)z^{-n}$
$H_{3}(z) = \sum h(4n+3)z^{-n}$

To do this in matlab i just did\

 h0=h(1:4:end);
 h1=h(2:4:end);
 h2=h(3:4:end);
 h3=h(4:4:end);

This doesn't seem to produce the right results.... If anyone can clarify how to do polyphase decomposition in mat lab or tell me what I am doing wrong I would appreciate it

EDIT:

I have a signal x[n] that is a FDM (frequency Division Multiplex) signal made up of four real valued baseband signals interpolated by a value of of 4 and modulated using the complex carrier $e^{j \frac{2*pi*k}{4}n}$  thefour base band signals are added to make the FDM signal x[n].

My goal is to estimate each of the 4 base band signals using the polyphase components of the lowpass filter.  The matlab code I am using is:

h = firls(100, [0 ,2 ,.25 1], [ 1 1 0 0],[.001 60]);
h0=h(1:4:end);
h1=h(2:4:end);
h2=h(3:4:end);
h3=h(4:4:end);

%first baseband signal
y0 = conv(x,h0);
y0 = downsample(y0, 4);
y0 = abs(fft(y0));

Here are the 4 estimated baseband signals I get compared to the actual

The polyphase components don't seem to be changing anything

I also tried setting up the filters like this

% Lowpass Filter
h0 = firls(100, [0 .2 .25 1], [1  1 0  0],[.001 60]);
h1 = firls(100, [0 .2 .25 1], [1  1 0  0],[.001 60]);
h2 = firls(100, [0 .2 .25 1], [1  1 0  0],[.001 60]);
h3 = firls(100, [0 .2 .25 1], [1  1 0  0],[.001 60]);
figure(1)
impz(h);
figure(2);
freqz(h);

% PolyPhase Componets
h0(2:4:end)=0;
h0(3:4:end)=0;
h0(4:4:end)=0;

h1(1:4:end)=0;
h1(3:4:end)=0;
h1(4:4:end)=0;

h2(1:4:end)=0;
h2(2:4:end)=0;
h2(4:4:end)=0;

h3(1:4:end)=0;
h3(2:4:end)=0;
h3(3:4:end)=0;

and ended up with these plots

still seems to be something wrong

(NOTE: Orange is actual signals.  Blue is my estimate.)


What is the relationship between group delay and propagation delay for a filter? I was trying to design a filter for a system and I was able to obtain the group delay plot in matlab using grpdelay. The group delay output is in samples. However, I was not able to relate that data to the propagation delay of the signal through the filter. Is there a formula between these two quantities?
Also, if I were to change the magnitude of the input signal in a step (keeping the frequency constant), how long would this change take to propagate to the output? Would that depend on the group delay or phase delay or would it depend on the filter order?
For example, if the filter cutoff is 5000Hz and I have unity gain at 500Hz and Group delay of 15 samples and Phase delay of pi radians, then if I send a 500Hz signal through the filter, after how many samples would I get it at the output? Assuming sampling frequency is 10000Hz.

The quantity you might be looking for is signal-front delay, which is the delay of the beginning of a signal passing through a linear system. It is simply the largest value $\tau_{sf}$ (in samples) for which
$$h[n]=0,\qquad n<\tau_{sf}\tag{1}$$
is satisfied, where $h[n]$ is the system's impulse response.
If an input, or change in the input, starts at $n_0$, then the system's response to it begins at $n_0+\tau_{sf}$.
If the impulse response has a long period of pre-ringing, the signal-front delay is generally much smaller than the actual propagation delay. In that case, definition $(1)$ can be modified as
$$\big|h[n]\big|<\delta,\qquad n<\tau_{sf}\tag{2}$$
where $\delta$ is a positive constant which is small compared to the maximum value of $h[n]$.
Depending on the nature of the impulse response, another meaningful definition of signal delay introduced by a linear system could be the center of gravity of the impulse response:
$$\tau_{gr}=\frac{\displaystyle\sum_n nh[n]}{\displaystyle\sum_n h[n]}=\frac{jH'(0)}{H(0)}\tag{3}$$
where $H(\omega)$ is the system's frequency response, and $H'(\omega)$ is its derivative w.r.t. to $\omega$. Note that $\tau_{gr}$ is simply the group delay evaluated at $\omega=0$.
Using the center of gravity as an estimate for signal delay is only useful for systems with a lowpass characteristic. For other systems (such as highpass or bandpass filters), the denominator of $(3)$ could be zero or very close to zero. For that reason, for general systems one would exchange $h[n]$ in $(3)$ for either its magnitude or its square.
Note that unlike phase delay and group delay, neither signal-front delay (including its modified form $(2)$) nor the center of gravity are functions of frequency.


Estimate Arrival Time of a Signal I'm looking for a Python package or method to accurately estimate the arrival time of time signals.

A single time signal is stored inside a 1-dimensional numpy array, so every component of the array represents the value of the signal at a given timestep.

Here is an example of the signals shape:

Setting a threeshold is not always efficient, it depends on the signal behaviour and amplitude.

I would like to find a gradient based method to identify the starting timestep.
Do you have any relevant experience or suggestion?

A gradient-based method is essentially a high-pass filter, if you think about it. ("gradient" being "derivative" for 1D signals, and the derivative being small for non-changing signals, and very high for high-frequency signals). And that's exactly how I implement it:

Signal in
high-pass filter
abs()
find first above threshold
subtract delay of filter¹

You can do so by designing a high-pass filter using scipy.signal.fir_filter_design and apply it using scipy.signal.filter.

Another popular method, which is mainly interesting if you want to do this on a live stream of samples, is doing the same within GNU Radio, where you can take the filter taps from the filter design, and put them into a readily available FIR filter. I kind of covered that in a blog post once, and you might want to read the first chapters of the GNU Radio Guided Tutorials, too.

¹: this requires your filter to have constant group delay, i.e. be linear-phase.


Periodicity of peaks within a signal I have a signal, in which i want to detect the period of the peaks, that exist within this signal: for example : After each 5mm or (5 seconds depending on the sampling rate and speed of my motor) i have a peak in my motor current.
I though of using the FFT, but then i had some problems with the sampling rate which lead me to look for the ununiform fft, but i wasnt very successfull.
I also tried using the Savitzky golay filter to smooth the data, since i read that it actually helps finding peaks within the signal, but not the period.
I tried appliying the FFT to the smoothed data, but I wasn't able to find any frequence.
So i was wondering is there is an algorithm that detects the periodicity of peaks within a signal.
 this is the Normal Signal plot .

 The signal in display(20mm/s) has been smoothed using the Savitzky golay filter, with a window of 51.

Before applying the algorithm below, you should first apply a DC blocking filter. That should remove the DC offset (and partially the step at the beginning), which will allow more usual periodicity techniques to be used.

You can try using the average square difference function. There is a nice write up on this answer and a way to implement it in an efficient manner here.

Here's an example plot. The $\color{red}{\bf red}$ circle shows the selected minimum. Note how the "octave" problem shows up: if the signal is periodic with period $P$, it's also periodic with period $N P$ where $N$ is an integer.

R Code Below

#30511

pulse <- 16 - seq(-4,4)^2
len <- length(pulse)

T <- 1000
P <- 49

ix <- seq(1,T - len, P)

data <- rep(0,T)

for (this_one in ix)
{
  data[this_one + seq(0,len-1)] <- pulse
}

data <- data #+ rnorm(T,0,1)

plot(data, type="l")

asdf <- function(data, Pmax)
{
  N <- length(data)
  asdf <- rep(0,Pmax)
  for (p in seq(1,Pmax))
  {
    ix <- seq(p+1,N)
    asdf[p] <- sum((data[1:length(ix)] - data[ix])^2)
  }

  return(asdf)
}

calc <- asdf(data, 5*P)

par(mfrow=c(2,1))
plot(data, type="l")
title('Repeating pulses')

plot(1:(5*P), calc, type="l")
period <- which.min(calc)
points(P, calc[P], col="grey", lwd = 12)
points(period, calc[period], col="red", lwd = 5)
title('ASDF')


MATLAB: Applying Inverse Fast Fourier Transform to an array of real numbers I'm trying to understand the difference between applying the Discrete Cosine Transform (DCT) and the Inverse Discrete Fourier Transform (IDFT) to the log Mel-filterbank energies as explained in the answer here. So, I tried the following pretty simple example in MATLAB:

x1 = [1. 2. 3.];
X1 = fft(x1);
x2 = ifft(abs(X1));
x3 = ifft(x1);

The classes of the results are:

Why do x2 and x3 have two different classes?

Update:

Here is the file 'myVoice.wav' in the following test:

[signal, fs] = audioread('myVoice.wav');
segmentLength = 240;
x = signal(20490 : 20490 + segmentLength - 1);
x = x .* hamming(segmentLength);
res = ifft(abs(fft(x)) .^ 2);

Since $\vert X_1 \vert$ is a real even function, its FFT (IFFT) is a real function. This is a basic property of Fourier transforms.

The power spectrum can be real but not even, in which case the ifft will give complex coefficients. In the example you linked, they start with a non-even real signal and then take the FFT. This will give complex values coefficients for the FFT. The answer then says that

  The problem is that we have only our log-energies, no phase
  information, so after applying the ifft we get complex-valued
  coefficients

So, after taking the FFT, the phase is discarded, and this loss of information makes it impossible to get back to the original real signal. (The Fourier transform is bijective). Once the phase is discarded you have a real but non-even function, the IFFT of which will yield complex coefficients.


In filter design, why isn't it possible to have a frequency response with phase 0? A linear phase filter will cause a delay in the output signal of the filter.

If I'm not mistaken, a filter with phase equal 0 would cause no delays in the output signal.

Can someone explain me why isn't it possibly to generate a 0 phase response filter? I mean, theoretically...

I'm studying for a test and I can't find the answer in textbooks.

A zero-phase frequency response is real-valued, and, by a property of the Fourier transform, the corresponding impulse response is symmetric. Since it is symmetric (i.e. it is non-zero for $t<0$) it can't be causal, and consequently, it can't be implemented (without adding delay, which defeats the purpose).


How do I get this Fourier transform? In the image, the signal is a pseudorandom sequence with values oscillating between -1 and +1. It is periodic (in this case with period 0.001). I don't see how the sequence yields a fourier transform whose real and imaginary parts are given on the right side.
Any help on how to compute this would be really appreciated.

You did not specify what sampling frequency you used to compute Fourier Transform. Software like MATLAB/GNU octave are only simulating an analog signal with high enough sampling rate. As a very general rule of thumb, if you know that time domain signal has a high periodic component of 1KHz, a 10x sampling frequency will give comprehensive results.
The pseudo-random waveform you showed on the left is a real and periodic signal and so we see real part of the FT is even-symmetric and imaginary part is odd-symmetric essentially meaning you can drop the left half plane from FT and focus on frequency values greater than zero. If you choose a sampling frequency correctly, magnitude response of the FT should show you peaks on frequencies which have the major share in your time domain signal e.g. 1KHz and 3KHz cosine signals may makeup (of course along with many smaller contributions from other cosine signals) the waveform. The real and imaginary parts are useful to compute the phase plot from FT essentially telling you how much delay should be added to cosine signals in order to make an approximation of your time-domain waveform.


Interpolation of magnitude of discrete Fourier transform (DFT) For example for peak frequency finding, it seems valid to use band-limited interpolation methods on the complex DFT bins, or separately on their real and imaginary parts and to calculate the magnitudes or squared magnitudes of the results. But how about band-limited interpolation of the magnitudes of the bins (I don't think that's valid), or their squared magnitudes (maybe valid)? By valid I mean that perfectly interpolated values should be equal to those found by calculating them from a larger DFT of a zero-padded version of the time-domain signal.

The first approach guarantees a non-negative result, unlike the others if the interpolation is not perfect, see this question about nonnegative or positive band-limited interpolation.

First a demonstration that the squares of both
$$\begin{align}&[\dots, 0, 0, 1,\hphantom{-}1, 0, 0, \dots] \text{ and}\\
&[\dots, 0, 0, 1, -1, 0, 0, \dots]\end{align}$$

equal

$$[\dots, 0, 0, 1, \hphantom{-}1, 0, 0, \dots]\hphantom{\text{ and}}$$

but the squares of their sinc interpolations differ (Fig. 1):

Figure 1. Squares of sinc interpolations of $[1, 1]$ (blue) and $[1, -1]$ (red).

This demonstrates that it is not in general possible to recover the square of a band-limited signal from its uniform samples taken at the critical sampling frequency of the band-limited signal.

Let's test different interpolation approaches in Octave. The gold standard of band-limited interpolation is DFT–zero-pad–DFT:

>> format free
>> x = [1 2 3];
>> y = [1 2 3  0 0 0];
>> abs(fft(x))
ans =

 6 1.73205 1.73205

>> abs(fft(y))
ans =

 6 4.3589 1.73205 2 1.73205 4.3589

The last set of numbers are the magnitude values calculated from perfectly interpolated frequency domain bins. Let's try interpolating magnitude instead:

>> fft(fftshift(horzcat([0 0], fftshift(ifft(abs(fft(x)))), [0])))
ans =

 6 4.57735 1.73205 0.309401 1.73205 4.57735

Seems quite a bit off. Now let's try interpolating squared magnitude:

>> fft(fftshift(horzcat([0 0], fftshift(ifft(abs(fft(x)).^2)), [0])))
ans =

 36 25 3 -8 3 25

>> sqrt(ans)
ans =

 (6,0) (5,0) (1.73205,0) (0,2.82843) (1.73205,0) (5,0)

It's not only off, but also sqrt() returned a complex number for a negative interpolated value. So is the only valid way to interpolate the bin values?

Let's give this one more chance, by trying to interpolate frequency domain data upsampled by a factor of 2. Because of the even length transform, this requires duplicating the "Nyquist sample", so apologies if the code is getting hard to read.

>> z = [1 2 3  0 0 0  0 0 0  0 0 0];
>> abs(fft(z))
ans =

 6 5.55485 4.3589 2.82843 1.73205 1.77302 2 1.77302 1.73205 2.82843 4.3589 5.55485

The above is what we want. Let's try interpolating 2x upsampled magnitude:

>> fftshift(horzcat([0 0 0], fftshift(ifft(abs(fft(y)))), [0 0 0]))
ans =

 3.36365 1.10447 0.318175 0 0 0 0 0 0 -0.208949 0.318175 1.10447

>> ans(4) = ans(length(ans)-2)
ans =

 3.36365 1.10447 0.318175 -0.208949 0 0 0 0 0 -0.208949 0.318175 1.10447

>> fft(ans)
ans =

 5.79105 5.59483 4.56785 2.7273 1.5231 1.76882
 2.20895 1.76882 1.5231 2.7273 4.56785 5.59483

That is still off. Let's try interpolating 2x upsampled squared magnitude:

>> fftshift(horzcat([0 0 0], fftshift(ifft(abs(fft(y)).^2)), [0 0 0]))
ans =

 14 8 3 0 0 0 0 0 0 1.18424e-15 3 8

>> ans(4) = ans(length(ans)-2)
ans =

 14 8 3 1.18424e-15 0 0 0 0 0 1.18424e-15 3 8

>> sqrt(fft(ans))
ans =

 6 5.55485 4.3589 2.82843 1.73205 1.77302 2 1.77302 1.73205 2.82843 4.3589 5.55485

Now it works perfect! The take home message is to upsample (time domain zero pad) at least by a factor of two before trying to interpolate in the frequency domain, and to interpolate squared magnitude rather than magnitude. It works because taking the squared magnitude is the same as multiplying each bin value by its complex conjugate. Complex conjugation preserves the bandwidth of the bandlimited function represented by the data, so the multiplication doubles the "time domain bandwidth" because it is equivalent to time domain convolution. Note that when choosing the interpolation method the 2x upsampled squared magnitude is still critically sampled, so further oversampling should make accurate interpolation much easier.

I about answered my own question, but will accept further insight as an answer!

P.S. Just found out there is also interpft, which does the interpolation with less syntax.

Utilizing additional information

Interpolation becomes easier or even exact with additional information about the data, for example that it is a critically sampled squared time-shifted sinc. In that case, given the two samples $\alpha$ just before and $\gamma$ just after the largest sample, the time $-1<d<1$ of the peak can be calculated by:

$$d = \begin{cases}0&\text{if }\alpha = \gamma,\\
\frac{1 - \sqrt{1 - p^2}}{p},\, p = \frac{γ - α}{α + γ}&\text{otherwise,}\end{cases}$$

with $d = 0$ meaning that the sinc is shifted to exactly the time of the largest sample. The same interpolation can be done at half sample rate, which is the critical sampling frequency of the underlying function the magnitude of which equals that of a time-shifted sinc, with the two successive largest-valued samples $\alpha$ and $\beta$ of the square of the absolute value of the underlying function:

$$0 \le d \le 1\\
d = \begin{cases}0&\text{if }\beta = 0,\\
1&\text{if }\alpha = 0,\\
\frac{1 + p - \sqrt{1 - p^2}}{2p},\, p = \frac{\beta - α}{α + \beta}&\text{otherwise,}\end{cases}$$

The formulas are not affected by amplitude scaling of the data. For frequency estimation you'd rarely have a purely real time-shifted sinc, but if you do, there are exact interpolation formulas for it.


Large sample rate or large number of cycles when estimating frequency? I wish to determine the frequency, amplitude and phase of a sine wave with some noise and I can only have a fixed number of samples. Should I go for many cycles and a small sample rate or a few cycles and a large sample rate?

I know that if I choose to use a DFT to determine the frequency then the spectral resolution will be the sample rate divided by the number of points. I must avoid aliasing so this suggests I should choose a small sample rate giving for example 3 or 4 points per cycle and capture many cycles of data. It may be a separate issue but the actual frequency will probably not correspond to a line in the spectrum so there will still be the, I think, unresolved issue of finding the frequency amplitude and phase from the small cluster of points around the maximum in the spectrum.

I could also choose an alternative method where I curve fit a sine wave to the data. This should give me best estimates in the least squares sense. However, for this method it is not clear which is best many cycles with a small sample rate or fewer with a large sample rate.

Is this a known problem with a known solution? Thanks for any input.

For frequency estimation of a pure sine wave, the difference between 2 slightly different frequency sine waves increases over time.  So a longer time base (as long as the sample rate is well separated from the Nyquist rate, or any integer sub-multiple thereof) would provide a bigger difference in a fixed number of samples, thus making it more likely that two different frequency signals could be separated from the effects of noise on any one signal.  This bigger separation difference is a close equivalent to better frequency estimation accuracy.  Note that under-sampling is possible if the S/N ratio is good enough, and as long as you know, a-priori, the correct frequency range of your signal and stay away from integer multiples of Fs/2 being near to Festimated.

A longer time base also corresponds to a narrower frequency bin spacing of a DFT result, thus decreasing the percentage errors in any interpolation between bins that may be performed or required.  But any frequency is a line in the spectrum after (high-quality/Sinc) interpolation.

Once you estimate the frequency, you can use two or more points (depending on the S/N ratio) to least squares fit a sine wave equation in two unknowns to estimate amplitude and phase.  Or windowed-Sinc interpolation, followed by successive approximation is another method to estimate amplitude.  Quadrature modulation of the signal against the estimated frequency can also be used to estimate phase (most accurately when referenced to the point or time at the center of the data window).


When applying a bandpass filter designed in Matlab's DSP system toolbox, is it necessary to do any normalization? I'm applying the next bandpass filter to the data of an impulse response coming from .wav; it fades data below 10 Hz and above 400 Hz:

d = fdesign.bandpass('N,Fst1,Fp1,Fp2,Fst2,C',...
                117,9.99,10,400,400.01,10000);
d.Stopband1Constrained = true;
d.Astop1 = 60;
d.Stopband2Constrained = true;
d.Astop2 = 60;

Hd = design(d,'equiripple');    % Passband type
h = fvtool(Hd);    % Filter design plot

Whose plot is:

Original signal vs. filtered signal looks like this:

Resulting signal is pretty much different from original, so I took a look into frequency domain:

As you can see, magnitude order is notoriously smaller in the bottom plot than the top one in this second plot. Considering the first picture, at least for around 200 Hz resulting amplitude might be almost equal to the original one but it doesn't happen (regardless of the amplitude increase after 300 Hz).

Is something missing in my design that causes this differences? Do I have to do any normalization to get an amplitude after filtering similar to the original one? Using 117 as filter's order is the best option to minimize fading.
Thank you in advance for you help.

Note that you are cutting a lot of signal frequencies in the transition from 10Hz to 200Hz (look the attenuation in the frequency response of the filter), consider to use lowpass filter or to change filter prototype to someone with faster transition.
By cut definition each frequency that is attenuated more than -3dB is cut.


Simple discrete time transformation $x[(n-1)^2]$ If I want to do the following transform
$$x[(n-1)^2]$$

I have thought about substituting each value of $n$ in $x[n]$ according to $x[(n-1)^2]$

i.e at $n=0,    x[0]=1$  to put it in $x[(n-1)^2]$

$x[(0-1)^2] = 1$ and so on...

Am I right ? and if I am at $n=-4$ it will be replaced at $n=25$


piecewise-quadratic vs piecewise-cubic vs higher order polynomial interpolation? There is a question available on DSP SE that mentions types of interpolation used for signal reconstruction but there isn't any mention about the difference between piecewise-quadratic,piecewise-cubic and higher order polynomial interpolation
So how we can differentiate between them ?
Link of question is given below:
Types of interpolation used for reconstruction in DSP?

The mathematically correct way of signal reconstruction would be sinc interpolation https://en.wikipedia.org/wiki/Whittaker%E2%80%93Shannon_interpolation_formula . However that's not very practical so most real world implementation use an approximation. But it's still based on the same principal: it is essentially a low pass filter and NOT a piecewise polynomial interpolation.
That's how most of real world signal reconstruction is done.
Other interpolation methods may be useful in certain special cases for specific applications. To differentiate between these methods, you would have to formulate what application specific metrics distinguish between a "good" and a "bad" interpolation. Then you can analyze different piecewise interpolation techniques with "typical" signals using the defined metrics and see how the performance is impacted and what is optimal.


feature detectors and descriptors comparison There are several kinds of detectors and descriptors, like SIFT, SURF, FAST.
I wonder are they all eligible for real-time applications? Which is the best or better?

And furthermore, is Harris-Laplacian dectector still useful when we already have the above three? Is it better than them?

I think SURF is better than all , as its clear from its name , Speed Up robust Feature Detector , but its not free to use , its a patent , so you can't use it for your commercial uses , you can only use it for your studies and knowledge purpose , but i think its also depend on your app that algo to use , but surf is better with speed


Cross-Correlation : Comparing two signals or finding the location of a target I read across places that cross-correlation of two signal (a radar sending a signal and receiving it back) one can detect the location (distance) of a target from a point and that cross-correlation can be used for comparing signals.

I solved a question to find cross-correlation of two finite length sequences where:

x(l)={1,2,1,1}
y(l)={1,1,2,1}

The cross-correlation sequence came as:

{1,4,6,6,5,2,1}

The answer was also correct as per my verification.

But I am not able to understand how using the above cross-correlation sequence value I can state it to be used for -

  a)location of a target? b)comparing signals x and y?

Your example would be more illuminating for you if you did this:

x(l)={1,2,1,1}

y(l)={0,0,0,0,0,0,0,1,2,1,1}

The peak of the crosscorrelation would then tell you how far into y(1) the contents of x(1) are.  If x(1) represents a radar pulse, then y(1) represents the received reflection.  The peak of the cross correlation then tells you how long it took for the reflection to return.  This tells you how far away the obect that reflected your signal - you then know where that object is.

Cross correlating also tells you how similar two signals are.  The higher the peak, the more the two signals are alike.

Try these to see what I mean:

These have a relatively low peak

x(l)={1,2,0,0}

y(l)={0,1,1,1}

compared to this next set which has a relatively high peak:

x(l)={1,2,0,0}

y(l)={1,2,0,0}


Harmonics of square wave I have a square wave, with 10% duty cycle. It is 90ms LOW, 10ms HIGH. Sampling frequency is 1kHz.

In spectral domain, after following this example, I get fundamental frequency (10Hz), as well as its harmonics. Also, the whole spectrum is enveloped with harmonics if 100Hz (which comes from pulse duration 10ms)

That is simply the case of spectral decomposition working properly for square wave, because of Fourier transform definition, right? I also read similar question but still a bit confused

Here is one way to think about this:

Start with the rectangular pulse of 10ms length. DFT is the sinc() function with a "period" of 100Hz
Now repeat periodically with 100ms intervals. Mathematically that's a convolution with a 100ms pulse train.
The DFT of the pulse train is also a pulse train in the frequency domain with 10 Hz spacing
Convolution with in the time domain is multiplication in the frequency domain. Hence we need to multiply the 100 Hz sinc() with a 10 Hz pulse train. That's the basic shape you have: a line spectrum with 10 Hz line spacing with the overall contour of a 100Hz sinc()
Sample at 1kHz. The sinc() function decays fairly slowly with frequency: around the Nyquist frequency you are only at -25 dB or so. Hence you will get non-trivial aliasing. However, since your sample rate is an integer multiple of the line spacing, the aliasing products will fall on top of the existing lines and will be hard to see. It'll just be a small deviation from the sinc() contour.

The "fuzz" between the lines appears to be noise. Your time domain wave form looks noisy and not particularly clean. It could also be caused by the way you calculate the spectrum (windowing, periodicity misalignment, FFT size, etc.)


How to prevent "zipping" effect on a modulated fractional audio delay line? (Flanger) I am implementing a Flanger using a fractional delay line.  I am modulating the length of the delay line using a sin function.  The delay line already uses linear interpolation to compute the delay output for the fractional sample.

When the delay length is decreasing (and pitch is increasing), the modulated length sounds just about perfect, with no artifacts.  However, when the delay length is increasing (and pitch is decreasing), I am getting a kind of "zipper" effect, with intermittent clicks, that one would expect from a non-interpolated delay line, as if sharp edges are being produced.

My guess is that I need to do additional interpolation beyond the simple fractional interpolation, but I have some questions.

First, why would I only hear artifacts during the period where my modulated delay length is decreasing?

Second, what is the best method to reduce these artifacts?  My first instinct is an all-pass filter, but is it typical to use an all-pass filter as well as an interpolation over a fractional delay length?  I just want to make sure I'm getting it right, and not wasting unnecessary cpu cycles :).

Some code will help.  In this code, config holds all delay line info, and delaySamples is the current working length of the delay line.  There are two read pointers, two write pointers, and two feedback values - one for each channel.  Thanks for taking a look!

void delayconfig_setDelayTime(DelayConfig *config, float delay) {
    config->delayTime = delay > 0 ? (delay <= 1 ? delay : 1) : 0;
    if (config->delayTime < 0.0001) config->delayTime = 0.0001;
    config->delaySamples = config->delayTime*SAMPLE_RATE;
    int channel;
    float *rp, *wp;
    pthread_mutex_lock( &delayMutex);
    for (channel = 0; channel < 2; channel++) {
        rp = &(config->rp[channel]);
        wp = &(config->wp[channel]);
        *rp = *wp - config->delaySamples;
        while (*rp < config->delaySamples) {
            *rp += config->delaySamples;
        }
     }
     pthread_mutex_unlock( &delayMutex);
}

float interp(float *delayBuffer, float delaySamples, float position) {
    position = fmod(position, delaySamples);
    int prev = floorf(position);
    int next = fmod(ceilf(position), delaySamples);
    float frac = position - prev;
    return delayBuffer[prev] + frac*(delayBuffer[next] - delayBuffer[prev]);
}

void delay_process(void *p, float **buffers, int size) {
    DelayConfig *config = (DelayConfig *)p;
    float out;
    int channel, samp;
    float *wp, *rp;
    for (channel = 0; channel < 2; channel++) {
        rp = &(config->rp[channel]);
        wp = &(config->wp[channel]);
        for (samp = 0; samp < size; samp++) {
            // here is where I'm changing my delay length.  The values are tentative while I'm figuring this out, and aren't really improtant.
            // The important thing is thing is that it's a sin function
            delayconfig_setDelayTime(config, 0.0061f + 0.006f*sin(M_PI*config->count++*INV_SAMPLE_RATE));
            pthread_mutex_lock(&delayMutex);
            if (*rp >= config->delaySamples) (*rp) -= config->delaySamples;
            if (*wp >= config->delaySamples) (*wp) -= config->delaySamples;
            float interpolated = interp(config->delayBuffer[channel], config->delaySamples, (*rp)++);
            int wpi = floorf((*wp)++);
            pthread_mutex_unlock(&delayMutex);
            out = interpolated*config->wet + buffers[channel][samp]*(1 - config->wet);
            if (out > 1) out = 1;
            config->delayBuffer[channel][wpi] = buffers[channel][samp] + out*config->feedback[channel];
            buffers[channel][samp] = out;
        }
    }
}

It's a bug.

M_PI*config->count++

You are processing 2 channels, but incrementing your modulation in the inner loop. So at the end of processing 1 channel, you continue to increment the modulation as you process the second channel.

c1 [ 1,2,3,4]         [9,10,11,12]...
c2           [5,6,7,8]

The modulation has a discontinuity, i.e in my simplifies example - c1 modulation jumps from '4' to '9'. resulting in zipper noise. Hope that makes sense.


Conceptually, how does real time convolution reverb work? I've got an impulse response of a hall, it's 10 seconds long.

Now say I want to apply this reverb to a wav file of some singing. I would perform a convolution of the data in two wav files and that the result would be a third lot of data, this time with the reverb. Great.

But how would this actually work in a real time setting, in a program like GarageBand, using a live input like an electric guitar?

Garageband can't convolve the full guitar track because it doesn't yet exist (because it's being played live!) One way that I thought of, was if the program performed the convolution on say, the 256 most recent samples of guitar input.

But this raises 2 issues:

Would the computer be able to perform so many convolutions every second?
And surely this introduces some sort of delay while waiting for the calculations to be processed?

If anyone could help explain, that would be great.

It comes down to latency vs. complexity. If your filter is 10 seconds long, you need to store the audio data of the last ten seconds and then you are able to calculate the current output audio sample with a latency of basically zero (ignoring the time required for calculations here) simply by doing:

$$y[0] = \sum_{k=0}^{l} x[-k] \dot h[k]$$

where $l$ is the length of the impulse response $h$ and $x[0]$ is always assumed to be the current sample. While this gives you basically zero lag, it boils down to the performance of naive convolution which can easily be too much especially if the impulse response is very long.

Now you could use techniques such as overlap-save which already perform convolution of very long signals by splitting the signal up and performing the convolution in chunks. However, to be able to calculate the output you first need to capture all the input samples for one chunk so this defines your latency (again ignoring the calculation time). The efficiency of these algorithms is usually in between the direct convolution method and the pure FFT method and the amount of calculations per sample will be less for longer chunks at the cost of latency.

The problem with overlap-save is that the chunk needs to be at least the size of the impulse response so for a 10 second response, the latency will be at least 10 seconds, hardly real-time. For this case so-called partitioned convolution algorithms have been developed that split both, the impulse response and the incoming signal. I've found a very good paper describing these techniques here (Click me) so I will just link it. I can hardly describe it better than they do.

From what I heard, there should be no more than 3-5ms latency or a good musician will notice the delay.


Concepts Surrounding First Two Samples in Difference Equation I'm attempting to construct a difference equation in matlab for the function u(n) and find the auto-correlation matrix R.

I'm following the technique as is defined by my text below:

The system I'm trying to model is also defined by the text below. My a1 and a2 are different.

I have it defined in Matlab as is shown below and included results for my auto-correlation matrix. Let me clarify that I know that technically R is wrong because I'm not computing the expected value. I hadn't figured out how to do that yet and am attempting to tackle what I perceive to be the larger problem first.:

a1 = 0.1;
a2 = -0.72;

m_max = 101; % Sample Length

v = randn(1,m_max); % Create m samples of Gaussian white noise
                % with 0 mean and unit variance

u = zeros(1,m_max); % Temp matrix for input vector

for n=3:m_max
    u(n) = v(n) - a1*u(n-1) - a2*u(n-2);
end

n = [1:1:m_max]; % Reconstruct n matrix to plot

u = transpose(u);
u_H = u';

R = u*u_H; % Corrleation Matrix

figure(1);
plot(v);
axis([1 m_max -5 5]);

This runs without any problem, but the issue is that my first two values are always zero. This makes sense in context of the code because I made the zero, but I know conceptually my correlation matrix R, as is shown by image 1 above,should not be 0 for r(0). Because my for loop thought must start from 3 (because u(n-2) is not defined for n<3) I'll always have 0's in my R matrix. I suspect that I would need to know initial values but those aren't provided to me by the problem. Am I thinking about this incorrectly? Any help would be much appreciated. Thanks!

This is just the transient effect of the filter due to zero initial conditions. The "solution" is to throw away the first two zero samples of $u[n]$ because they are not related to the input noise signal $v[n]$.


Why are the bass components of music more audible than other components? Consider a situation where I am listening to loud music in a room. If I move into another room and close the door (so music is now played in another room behind the closed doors), I will hear the bass component of the music much better than other components.

I would like to know why the bass voice is more audible (can be better heard through the closed door from another room) than other components in the music ?

The sound wave is a pressure wave in the air, a mechanical vibration. If you have two rooms separated by a wall or door, there is little or no air passage through which the sound can propagate.

Instead, the wall itself have to absorb and re-emit the mechanical vibration on the other side. Since doors and walls are massive objects, they are resonating on lower frequencies.

The high frequency sound cannot pass through the doors since it is able to make only small particles vibrating (or just surface of the door), hence the wave is not completely absorbed and rather reflected back.

The bass sound wavelenghts are about feet long (the range is quite large depending how you define "bass"), which roughly matches the size of such large objects like walls and doors.

Of course, waves with even larger wavelengths are unable to make objects of this size resonate and hence doors behave like a band-pass filter (as Dilip Sarwate noted).


What could cause a phase difference between a FPGA implemented FFT and numpy.fft.rfft We have a FPGA system which takes input from an ADC and calculates a FFT.  The system identifies the location of our frequencies of interest and sends the coefficients from those bins to my software.  In addition to the coefficients I am sent the time series data
While the hardware guys were working on the hardware I developed additional signal processing code.  To enable testing and development I wrote some code to synthesize the signals I expected to receive from the hardware and that works fine.  Not unexpectedly now that I have live data nothing works.
In the process of debuging the problem I am taking the real world time series data and using np.numpy.fft.rfft to look at the spectrum.  When I plot the spectrum I see our frequencies of interest in the correct FFT bins. I normalize the PSD and the Python's FFT and the hardware FFT match well.  The problem I see is the phases of the coefficients do not match.  The difference in phases don't look ordered.  (Just looking at the phase difference in a plot).
When I compute a FFT on the time series provided by the FPGA (this is the real world data) using numpy.fft.rfft I expect the coefficients at my frequencies of interest to have the same phases as the coefficients calculated by the FPGA FFT which is operating on the same time series.
Does anyone have an idea of what could cause FFT's on the same real world time series data to have different coefficient phases?
Thanks
Justin

Unless the sampling of the real signal is synchronized with the data, the sampled real data will often have different phases from synthesized data. The real signal will also not be noiseless.
Think of this as the difference between:
$$
x_{\tt real} [n] = A \sin (\omega n + \theta) + \epsilon[n]
$$
and
$$
x_{\tt synthetic} [n] = A \sin (\omega n )
$$
Here, $\theta$ is the phase offset due to the lack of synchronization and $\epsilon$ is the noise sequence.

Another thought: I'm wondering if it's a problem of offsets?  The manual entry for rfft says:

The truncated or zero-padded input, transformed along the axis
indicated by axis, or the last one if axis is not specified. If n is
even, the length of the transformed axis is (n/2)+1. If n is odd, the
length is (n+1)/2.

The numerical example at the bottom shows:

which indicates that while fft returns $N$ points, rfft returns $N-1$ points.  Could it be an indexing issue?


Attempting to use the trapezoidal rule to form a difference equation representing a circuit I have a differential equation that has been proven to be correct.

The transfer function obtained by Laplace domain analysis and Matlab freqs match up and all is well.

The problem is somewhere in my attempt to digitize the differential equation. I am required to do this using the trapezium rule to approximate the analog signal.

The goal is to get a digitized equation in the form of..

Vo^n+1 = alpha1(Vi^n+1) + alpha2(Vi^n) + alpha3(Vo^n)

Where each alpha value is made up of a different combination of T, R1, R2 and
C all numerical parameters.

I am following some guidelines set by my lecturer but if the approach seems wrong please let me know.

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

Current incomplete solution..

Then I am trying to convert that equation to the Z domain and therefore use the freqz tool in Matlab. Please see the working out below.

Here is the Matlab plot that shows what has gone wrong. The reason for the multiple lines are varying values of capacitance for C.

The blue lines represent the validated analogue freqs(laplace domain transfer function) frequency response.

The green lines show my attempt at using the freqz(z domain discrete transfer function) frequency response.

If it is of any further interest, the code to generate those plots can be found below as well.

Any help would be sincerely appreciated, clutching at straws and the occasional hair at the moment!

Cheers

MATLAB code:

%% Comparing the impulse amlitude response using two different methods
%
%  Basically checking if freqs matches freqz

clc; clear all;

% Circuit Parameters
R1 = 2.2e3;
C = 0.01e-6;
R2 = 10e3;

fs = 10000;   % Sampling frequency (Hz)
T = 1/fs;     % Timestep

% Coefficients
den = (T*R1 + T*R2 + 2*R1*R2*C);

N = 1000;

% Digital parameters
y = zeros(N,1);
y(1) = 1;
result = zeros(N,1);

% Analogue parameters
fa = 0:fs/N:fs-fs/N;
w = 2*pi*fa;
j =sqrt(-1);

lastValue = 0;

clf;
figure(1);

hold on;

for i=1:5

% Analogue Recalculations
b = 1/(R1*C);
a = (R1 + R2)/(R1*R2*C);
Ha = freqs(b,[1 a], w);

% Digital Recalculations
alpha1 = (T*R2)/den;
alpha2 = alpha1;
alpha3 = (2*R1*R2*C - T*R1 - T*R2)/den;

[Hd,fd] = freqz([alpha1 alpha1],[1 alpha3],N,fs);

for n=1:N-1

 currentValue = alpha1*y(n+1) + alpha2*y(n) - alpha3*lastValue;
 result(n) = currentValue;
 lastValue = currentValue;
end

C = 1.5*C;

plot(fa,abs(Ha),fd,abs(Hd)),xlabel('Frequency'),ylabel('H(s)');

end

hold off;

The Laplace transform of the original differential equation is
$$Cs V_o(s) = V_i(s) \left( \frac{1}{R_1} \right) - V_o(s) \left( \frac{1}{R_1} + \frac{1}{R_2} \right)$$
$$R_1 R_2 C s V_o(s) = R_2 V_i(s) - (R_1 + R_2) V_o(s).$$
In the Laplace domain, the ratio $V_o(s)/V_i(s)$ is the transfer function for the circuit,
$$H(s) = \frac{ V_o(s) }{ V_i(s) } = \frac{ R_2 }{ R_1 + R_2 + R_1 R_1 C s }
= \frac{K}{1 + R_p C s},$$
where $K = R_2 / (R_1 + R_2)$ is the voltage divider ratio for $R_1$ and $R_2$
in a series connection, and $R_p = R_1 R_2 / (R_1 + R_2)$ is the equivalent
resistant for $R_1$ and $R_2$ in a parallel connection.
We can obtain the $z$-domain transfer function for an equivalent digital filter
using the bilinear transformation,
$$s = \frac{2}{T} \left( \frac{z-1}{z+1} \right)$$.
In this transformation, the constant is set to $2/T$, which will produce a digital filter that is an approximation to the analog filter that is equivalent to approximating the convolution integral by the trapezoidal rule.  I assume this is what you are looking for.
By making this substitution, the transfer function for the digital filter is
$$H(z) = \frac{K}{1 + R_p C \left( \dfrac{2}{T} \right)
\left( \dfrac{z-1}{z+1} \right)}
= \dfrac{ KT (z+1) }{ T(z+1) + 2 R_p C (z-1) }$$.
This can be further simplified to
$$H(z) = \frac{V_o(z)}{V_i(z)} =
\frac{ \left( \dfrac{KT}{T + 2 R_p C} \right) \left( 1 - z^{-1} \right) }
{1 + \left( \dfrac{T - 2 R_p C}{T + 2 R_p C} \right) z^{-1}}.$$
Taking the inverse $z$-transform of this equation, we obtain a
difference equation for the digital filter,
$$Y(z) \left[ 1 + \left( \dfrac{T - 2 R_p C}{T + 2 R_p C} \right) z^{-1} \right]
= X(z) \left[ \left( \dfrac{KT}{T + 2 R_p C} \right) \left( 1 - z^{-1} \right) \right].$$
$$y[n] = \left( \dfrac{KT}{T + 2 R_p C} \right) \left( x[n] - x[n-1] \right)
- \left( \dfrac{T - 2 R_p C}{T + 2 R_p C} \right) y[n-1].$$
You should be able to use this equation to recursively compute $y[n]$.


How to calculate the Fourier transform of a mean filter in Matlab? In Matlab, how can I calculate the discrete-space Fourier transform of a mean which takes the average of 4 adjacent points, with this kernel

$$\begin{pmatrix}
      0 &1& 0\\
      1 &0& 1\\
      0 &1& 0\\
\end{pmatrix}$$

The 2D FIR 3x3 filter kernel you have provided is not finding the mean but the sum of samples. So you have to divide by 4 to get the mean. Furthermore it does not sum adjacent samples but diagonal ones.

Then if you create a 3x3 matrix h = [0,1,0 ; 1,0,1; 0,1,0] in matlab and you wish to perform a 2D-DFT on that matix h, then the following function call H = fft2(h) will return you the 3x3 discrete Fourier transform $H[k_1,k_2]$ (k1 rows, k2 columns) samples of the filter $h[n_1,n_2]$ as:

$$ H[k_1,k_2] = e^{ -j \frac{2\pi}{3} k_1} + e^{ -j \frac{2\pi}{3} k_2} + e^{ -j \frac{2\pi}{3} 2 k_1}e^{ -j \frac{2\pi}{3} k_2} + e^{ -j \frac{2\pi}{3} k_1}e^{ -j \frac{2\pi}{3} 2 k_2} $$ for $k_1,k_2 = 0,1,2$.

If you want to compute a larger $N \times N$ 2D-DFT on $h[n_1,n_2]$, then the following call H = fft2(h,N,N) will return you an $N \times N$ result as:
$$ H[k_1,k_2] = e^{ -j \frac{2\pi}{N} k_1} + e^{ -j \frac{2\pi}{N} k_2} + e^{ -j \frac{2\pi}{N} 2 k_1}e^{ -j \frac{2\pi}{N} k_2} + e^{ -j \frac{2\pi}{N} k_1}e^{ -j \frac{2\pi}{N} 2 k_2} $$ for $k_1,k_2 = 0,1,...,N-1$.


Stochastic Methods for Image Deconvolution Problem If we convolve an image with a point spread function and from the resulting image to find the input image, can we use any stochastic approaches? I feel like we will not be able to. A single image seems to me a deterministic quantity and I cannot think of any way to approach this deconvolution problem in a stochastic way. However, I am not sure and I want to know if there is a way. Any help is appreciated

Any Deconvolution method which takes into account the noise in the image is basically a stochastic approach.
Usually, the model for Deconvolution is:

So having the noise in there makes it a problem with stochastic properties.
Remark
If by stochastic you meant sampling from the Posterior Distribution then you may have a look at Stochastic Image Denoising by Sampling from the Posterior Distribution (Though it is not about Deconvolution but on Denoising).


Decoder implementation Im struggling to implement / build convolutional decoder .
I've done implementation of FEC encoder which it's as shown in the photo below:

So this is the encoder that I've implemented, and I want to do for it the decoder, from where could I start and how can I implement for this encoder the decoder? (decoder is working opposite to encoder and I know that but from where do I start and how to implement that decoder as what I did in my encoder photo)
thanks for any help/assistance.

It depends, but convolutional codes like these are usually decoded using maximum likelihood sequence estimation, typically through the Viterbi decoder.
The design methodology for any decoder doesn't start with an electronics implementation of the encoder (that is something that's highly specific to the speed and technology you want to use), but with a general design of the encoder, usually in shape of the polynomial that your circuitry implements.
Draw the state diagram to the encoder, then draw the trellis diagram to that, and then you gain something you can implement as viterbi decoder.
I can't see a case where someone would ever actually implement this encoder using discrete 74xx series logic gates like you do – I presume this is just a representation of the logic.
You can implement encoding logic very easily in software, running on a microcontroller; if you need 100s of megabits per second in throughput, you'd implement this in an FPGA, but then, this encoder wouldn't, from a code design point of view, be what you use to begin with, probably. So, sure, this is also something you can easily implement using logic gates.
Decoding is always more complex than encoding. Decoding convolutional codes is kind of inherently something that you'll need to test hypotheses while doing; this will be no fun / make no sense to implement in discrete logic components (gets way too large, and thus, too bug-prone very quickly), you'd do that in software, really, or on an FPGA.
Just like for your previous question: you can't understand how this works by implementing it in hardware first without reading any theory. So, wikipedia: convolutional codes, Trellis diagram, Viterbi decoder.


Signal correlation vs Pearson correlation between signals I have some data which are the sets of values uniformly distributed along some measure (simple example - along the time intervals).
I tested both the coefficients and found the difference, some times essential.
F.e. if both the data are noise-like, signal r is about 0.7, and Pearson r - about 0.25.

Which coefficient should be considered better and why?

I compute the signal correlation coefficient as
$$
r=\frac{\sum_ix_iy_i}{\sqrt{\sum_ix_i^2}\sqrt{\sum_iy_i^2}}
$$

and Pearson correlation coefficient as
$$
r=\frac{\sum_i(x_i-m_x)(y_i-m_y)}{\sqrt{\sum_i(x_i-m_x)^2}\sqrt{\sum_i(y_i-m_y)^2}}
$$

There are two fragments of pair of signals:

Fragment 1. Two signals with similar Signal- and Pearson- coefficients

Fragment 1 http://net2ftp.ru/node0/nf05@ngs.ru/Signals1.png

Fragment 2. Two noise-like signals: Signal r = 0.7, Pearson r = 0.2

Fragment 2 http://net2ftp.ru/node0/nf05@ngs.ru/Signals2.png

PS I do not need to performe cross-corelation, convolution etc.

Thank robert bristow-johnson for the right question.
In my case signal coefficient shows a some level of similarity between noise-like signals because of appreciable mean values. That is as you compare two electrical noises with some direct voltage: of course they are as more similar as high equal direct voltage they have.
Since I need to compare specifically shapes, I should use Pearson method only (which "normalizes" data itself in terms of mean values).


Auto-focus estimation? What kind of analysis might be used on sets of 2D monochrome image matrices or bitmaps to estimate or compare their relative quality of being "in focus"?  Is there an efficient way of computing a scalar metric on something like overall contrast, luminance sparsity, or other artifacts produced from out-of-focus blurring?

Images in focus have sharper edges, so applying an edge detector and measuring the energy of the output puts you on the right track.

A simple technique is to compute the sum of the laplacian over the images, possibly center-weighted if you deal with everyday photographs (people tend to put the object of interest in the center, so it's better to have the center in focus than the borders)

This criterion is for example what is used in many "Stack focusing" applications for macro photography.

There are other approaches used in de-hazing or tone-mapping... Search for "Local contrast measures".


RMS of absolute-valued sinusoid A signal $$x(t)=\left|\sin(2\pi(1000)t)\right|$$ ($\left|x\right|$ is the absolute value function) is fed to an ideal low-pass filter with cutoff frequency of $1500\,\mathrm{Hz}$ to produce an output $y(t)$, what is the RMS voltage of $y(t)?$

My initial belief is that it will be zero since the absolute value doubles the frequency of the signal and so the frequency of $y(t)$ will be $2000\,\mathrm{Hz}$ which is beyond the cut-off frequency of the ideal low pass filter.

However, upon reflecting on the Fourier transform of $y(t)$, I believe it will have a DC component with zero frequency which will give it an RMS voltage. Moreover, $y(t)$ is a fully rectified signal which, when passed through low pass filters in reality (i.e. capacitors) will yield a smooth DC signal.

How do I find the RMS voltage of the resulting $y(t)?$ Is it really nonzero? Will it be just the same as the RMS voltage of $x(t)?$

Your reasoning is correct. The signal $x(t)=|\sin(2\pi f_0t)|$ is periodic with period $T/2$, where $T=1/f_0$. Consequently, it has components at frequencies $2k/T$, $k=0,1,\ldots$, i.e., at $0$Hz, $2000$Hz, $4000$Hz, etc. Since the ideal low pass filter has a cut-off frequency of 1500Hz, all harmonics are completely suppressed, and only the DC component appears at the output.

I'm sure that you know how to compute the DC component of $x(t)$ (which, if positive, equals its RMS value).

Your concerns about actually implemented filters and capacitors etc. are irrelevant to the question, which explicitly states that the low pass filter  can be assumed to be ideal. But even a non-ideal low pass filter of appropriate order would suppress all harmonics sufficiently, such that they would be buried in noise.


Why does this constant-Q transform "lose octaves"? In this paper (Schörkhuber, Klapuri, Sontacchi), it's shown that the constant Q-transform evaluated "octave-wise" produces the following looking graph:

I don't understand why the octaves are "disappearing" in some time intervals. Isn't this then a lossy" process somehow?

In the paper it says that this is achieved because

  the atom hop sizes $H_k$ are set according to the highest frequency bin within each octave, that is, for each octave down $H_k$ is multiplied by 2

The lower frequency atoms are longer in time (or of a length roughly proportional to the period of the frequency, give or take 2X).  So less of these lower frequency atoms are needed per second to cover that time period.  And more of the shorter atoms (for shorter wavelength frequencies) are required for about the same coverage (or with roughly the same percentage overlap).

And to maintain a nearly constant Q, the length of a bandpass FIR filter ends up roughly proportional to some constant multiple of the wavelength of the center bandpass frequency of that filter.


Generalized/Power means in DSP Generalized/power means may be used to construct moving-average filters with different properties than regular one which is based on arithmetic mean. This observation seems to be trivial - even Wikipedia mentions that.
The problem is that power means are usually defined for sets of positive numbers - while in DSP we often deal with negative numbers too. Therefore, many of well-known inequalities and properties of power means may not be applied to sets with negative numbers.

This is of interest to me since I have found that in my partical area of application power means of higher order give better results than arithmetic mean. I have been able to prove some basic properties and I more or less understand where does the improvement come from. However, I would be interested in learning. However, everything I have found seems to focus on positive numbers. So my question is:
Are there any papers on using generalized means for moving average filtering or some other signal processing? Any papers on purely mathematical properties of generalized means for sets with negative values may be relevant too. Thanks!

Certain varieties of generalized means have been applied in signal processing, such as the Lehmer mean. This paper discusses several generalized means, mainly from a theoretical point of view.

Note, however, that all applications of generalized means deal with signals that are inherently non-negative. This is also the case in signal processing, where there are many applications with non-negative signals, such as frequency domain denoising, image processing, fuzzy logic, etc. The application of generalized means in image processing is discussed in this thesis.

I think that rather than looking for new definitions of generalized means for bipolar signals, it is probably more useful to define the signals properly such that generalized means can be applied.

A completely different alternative to generalized means might be general non-linear filters, such as Volterra filters, which don't have the restriction that the input signal must be non-negative.


What is this effect/filter in given images? Here are two images. What kind of filter/processing is done here?
It looks like an emboss but the bumps are smooth.

1) Original :

2) Processed :

Highlight blooming; the bright parts are bleeding into the shadows. It can be done with Photoshop's Diffuse Glow filter.


How to design a continuously time-variant digital filter? I have discrete time series containing signal with smoothly varied frequency over time (called a "sweep"). How can I design a discrete filter (low-pass or band-pass in my case) of a finite length with linearly varying cut-frequency over time and constant cut-slope?

EDIT: the signal is the sampled "trace" of the seismic source - a seismic vibrator, which sends the vibrations of the slowly varying frequency down the earth. The dependency of the frequency over time (the sweep) is known (let it be linear, $f(t)=f_1*(1-t)+f_2*t)$, but there is a problem that there might be another vibrators that operate on their own, and the task is to "band-guard" the trace of this vibrator avoiding the unwanted signals from other ones.

One approach would be to try to remove the frequency chirp from the observed data, thus translating all of the echoes to approximately baseband. I find this to be most straightforward by converting the observation to an analytic signal, then multiplying by a complex exponential whose instantaneous frequency is equal to $-1$ times the frequency chirp profile (while keeping its phase continuous). After dechirping the received data, you could then apply a lowpass filter to suppress any other sources that don't overlap in frequency with your chirp profile. If your followon analysis methods need to see the frequency ramp, you could reapply the chirp again by multiplying with another complex exponential.

The passband width of the lowpass filter defines how tightly around the transmitted tone that you reject other frequency components. The width of the passband would also need to be selected while taking the expected two-way propagation time of the transmitted signal; at time $t$, assuming a low-to-high frequency chirp, you might be transmitting frequency $f_t = f_c + \Delta f$, while the receiver is observing a delayed version of what you transmitted some time ago, for instance $f_r = f_c$. Your lowpass filter must have enough frequency coverage to cover the frequency slew of your chirp profile over the expected range of time delays. At the same time, though, you have an incentive to make the passband width as narrow as possible to reject other signal sources that are nearby your chirp profile in frequency, so as so often occurs in engineering, you have a tradeoff to examine.