How to calculate the impulse response of an RC circuit using time-domain method

Dynamic electric circuits involving linear time-invariant resistors, capacitors, and inductors are described by linear constant coefficient differential equations (LCCDE). Mathematical solution of such LCCDE requires some physical (electrical circuit theoretic) insight too. This is especially true for solving circuits under impulse functions (such as finding impulse responses). Note, that the dirac impulse $\delta(t)$ is a non-differentiable, non-continuous function; therefore it cannot be the subject of a continuous-differential equation. Nevertheless, a procedure as I present below, is developed to solve LCCDEs involving such improper impulse functions.
Given the circuit assuming zero initial conditions, derive its DE by applying a KCL (Kirchhoff's current law) at the output ($y(t)$) node :
$$ \frac{y(t)}{R} = C \frac{d}{dt} \left( x(t) - y(t) \right) $$
$$ y(t) = RC x'(t) - RCy'(t) $$
$$ \boxed{ y'(t) +\frac{1}{RC} y = x'(t) } \tag{1} $$
Impulse reponse $h(t)$ is the solution of Eq.(1) when the input is an impulse: $$ x(t) = \delta(t) \implies y(t) = h(t)  $$
Let's apply a two-stage procedure to find $h(t)$. First, modify RHS of Eq.(1) by replacing $x'(t)$ with $x(t)$ :
$$ y'(t) +\frac{1}{RC} y = x(t)  \tag{2}$$
and denote solution of Eq.(2) as $h_0(t)$ when $x(t) = \delta(t)$. Then, by the second stage, the overall impulse response $h(t)$ of Eq.(1) will be:
$$\boxed{ h(t) = h_0'(t) .} \tag{3}$$
Solution of first stage, Eq.(2):
From the characteristic equation : $s + \frac{1}{RC} = 0 \implies s = - \frac{1}{RC}$.
The homogeneous solution is :
$$y_h(t) = K e^{-t/RC} u(t) \tag{4}$$
The particular solution $y_p(t) = 0$, as the input is an impulse; i.e., $$\lim_{t \to \infty} x(t) = 0 \implies y_p(t) = 0$$
And the complete solution:
$$h_0(t) = y_h(t) + y_p(t) =  K e^{-t/RC} u(t) $$
To find the unknown $K$, we can use $y(0^+) = h_0(0^+)$ ; and from the circuit physics we find $y(0^+) = h_0(0^+) = 1$.
Then $K = 1$, and the total solution of Eq.(2) is
$$ \boxed{ h_0(t) = e^{-t/RC} u(t) } \tag{5}$$
$u(t)$ function is used to indicate that the output is causal; i.e., zero for $t < 0$. (a realizable physical circuit will be a causal system.)
Finally, the impulse response $h(t)$ is found from Eq.(3), by the second stage, as:
$$ h(t) =  h_0(t)' = \left( e^{-t/RC} u(t) \right)' $$
$$ \boxed{ h(t) = \delta(t) - \frac{1}{RC} e^{-t/RC} u(t) } \tag{6}$$
An intuitive physical explanation of the circuit behaviour under impulse excitation is as follows. Due to the physical capacitor's continuity of voltage  property, the initial impulse (physical spike) voltage applied at the input, appears right at the output as a result of KVL applied around the loop at $t=0$. This impulse voltage on the output resistor $R$ produces an impulse current $i(t)= \delta(t)/R$ at time $t=0$ looping clock-wise around the circuit. This impulse current passes through the capacitor and jumps its voltage from $V(0^-) = 0$ to $V(0^+) = 1/RC$. Then the input node is shorted as the impulse is gone at $t=0^+$. And the simple first-order capacitive discharge circuit is obtained. Note that the discharge current loops counter clock wise due to capacitor voltage polarity. Hence the output voltage becomes negative for $t>0^+$.


I've written articles like that, and it's incredibly hard -- and sometimes not the best exposition -- to get the graphs to match exactly with the text.  From the looks of things, the author was intentionally showing you idealized graphics rather than the actual output of code.

  My phase-difference signale_D[n] is flipped vertically from what the tutorial shows...

Eh, don't worry.  It's a control loop, there's always a stray minus sign in there someplace.  As long as you match it up with another stray minus sign so that it's stable, you're OK.

  I have output that is close to what the article shows, except that my cos_out[n] seems to be "leaning" to the left: cos_out

That's because your K_p is so high.  You can get a simple loop, fast settling, or a nice sinusoidal response -- not all three.  If you're using the author's numbers, or if you've selected gains to make your loop settle as fast as his idealized example (which was probably chosen more to make things fit in the page rather than because it's the best tuning) then what's happening is that strong signal at the sum frequency is bleeding into the NCO phase, and making the "lean".

Typically a PLL will be set up to take hundreds or thousands of cycles to settle, which implies a much smaller K_p (certainly less than 0.01, and almost arbitrarily small depending on what you're doing).  How fast a PLL settles depends on what you're doing, and what you're doing it with, so if your next question is "how fast should it settle" -- I can't answer.  But I would say try it for settling in around 100 cycles, 1000, and 10,000.  Keep in mind that the capture range (the frequency error between signal and NCO initial frequency over which it can actually acquire a lock) and its tracking performance (its ability to follow a varying input frequency) both go down with slower settling.


What is the input dimension for LMS update algorithm

The output of an FIR filter is
$$
y(n) = \sum_{i=0}^{N-1}w_i(n)x(n-i)
$$
Say the weight vector $\mathbf{w}(n)$ has a length of $N$, i.e.,
$$\mathbf{w}(n) = [w_0(n), w_1(n), \ldots, w_{N-1}(n)]^T$$
and the input vector $\mathbf{x}(n)$ should have the same length as $\mathbf{w}(n)$:
$$\mathbf{x}(n) = [x(n), x(n-1), \ldots, x(n-N+1)]^T$$
Therefore the output $y(n)=\mathbf{w}^T(n)\mathbf{x}(n)$ is a scalar, and the error signal $e(n)$ and disired signal $d(n)$ are also scalar. They represent the signals at discrete time $n$. You should update the filter weights sample by sample.


Implementation of a Variable Fractional Delay with Lagrange Interpolation using Farrow Structure

The task at hand is a pretty common one; most symbol timing synchronizers need to do exactly what you describe.

Same goes for many types of resamplers – they need to reconstruct varying fractional inter-sample "snapshots" from the input sample stream to construct a non-rationally related rate output.

So that's where I'd start looking.

In fact, the farrow structure idea takes you pretty far!

Imagine you could have a large bank of linear-phase filters. A linear-phase filter has a constant group delay.

By modulating said filters, you can change the phase slope and hence the group delay. Now, you just pick the filter that fits your needs.

In a polyphase filterbank, such as implemented by the GNU Radio "PFB Arbitrary Resampler" block, exactly that is done: You construct a whole bunch of filters, and for each sample, pick the one that's closest to what you need (or the two closest, and linearly interpolate between these).

Because it's used as a resampler, for every non-1.0 resampling ratio, the next sample will be choosen from a different filter branch, with the appropriate delay.

Now, since you, most of the time, just want to delay, you seldom need the "resampling" capability of that resampler. (Resampling is exactly what happens when you change the delay linearly)

So, instead, you'd have a single filter (if possible, use a filter that you need to have, anyway), and you'd modulate it.

Think about this in frequency domain: assume you had a linear, continuous Fourier transform that you apply to your signal and you multiply that with both the FT of your filter impulse response and $e^{j2\pi\tau f}$ before transforming back to time domain. To adjust delay, you adjust $\tau$. Tadaaah! Variable fractional delay.

Problem certainly is that you're in DSP, so you don't get that linear Fourier transform to implement linear convolution; you get the DFT, but multiplication in discrete frequency domain is equivalent to circular convolution in time domain.

But that's been solved already – fast convolution allows you to apply long filters with little effort (and low numerical problems – the FFT is exceptionally stable, usually even more stable than the direct convolution!). Now, if you go ahead and implement that e.g. via the overlap-save method, all you need to do is make sure that the phase of the complex exponential doesn't "jump" between FFT segments. You can then change (even gracefully!) the $\tau$ between individual convolutions.


Why is a feedback loop not represented by the least order transfer function?

Not really an answer, but let's work it out with $H(z) = h$ and $C(z) = \frac{K}{z-\alpha}$.:

$$
\begin{array}
\ L(z) &=& \frac{H(z)C(z)}{1+H(Z)C(z)}\\
&=& \frac{h \frac{K}{z-\alpha}}{ 1 + h \frac{K}{z-\alpha}}\\
&=& \frac{h K}{ z - \alpha + h K}\\
&=& \frac{h K}{ z - \alpha + h K}\ \ \frac{z - \beta}{ z - \beta}
\end{array}
$$

I added the last line to show what I mean by a pole-zero cancellation.  $\beta$ can be any value; if the pole and the zero at $\beta$ cancel out, the overall value of $\frac{z - \beta}{ z - \beta}$ is 1.

The trick will be to figure out why matlab inserts the extra $\frac{z - \alpha}{ z - \alpha}$ in there in the first plcae,


Are there any good "conceptual distances" for comparing two images?

Ok, I read the question posted in your link, so forget about MSE or SSIM.
What you can do is define a specific distance for the particular type of data you are using.
For example, if you are using video data, depending on what you are most interested in, you can define distance measure based on the average shot length, keyframe colour histogram, average motion in a shot, and so on.
If you are working on a database of face images, you can use e.g. PCA+eigenfaces to project each image on a N-dimensional space and compute an euclidean distance on that space; if you are working with e.g. flower pictures you can define a distance based on edge + colour information...

Here´s my previous answer:

It´s not really clear to me what you mean for robustness/transparency, to be honest.

If you want to check which one of I_2 and I_3 is more similar to I_1 a simple way is to use compute the mean square error, defined as the sum of the squares of all the pixels in the image, normalized for the total # of pixels


The pitch of a musical note is very often not the same thing as the frequency peak in an FFT magnitude computation.

This is due to all the overtones and harmonics within interesting musical sounds, as opposed to single pure (and boring sounding) sinewaves.  Please investigate pitch detection/estimation algorithms that do more than look for a bare FFT magnitude peak.  A little reading on human psychoacoustics, pitch perception, audiology, and the physics of musical instruments might also help.


Fourier series of $cos(\omega_0 t)$ in continuous time

They did not actually solve the integral, they just used a version of Euler's formula:

$$\cos(\omega_0t)=\frac12\left(e^{j\omega_0t}+e^{-j\omega_0t}\right)$$

which already is the Fourier series of $\cos(\omega_0t)$. Comparing $(1)$ to the formula for the Fourier series, you see that all coefficients are zero except for the ones with indices $k=1$ and $k=-1$, and they equal $\frac12$.


How do I show an impulse response leads to a zero-phase frequency response?

1) Define $d[n]$ as a sinusoidal

2) Apply angle addition formulas

3) Simplify

4) Interpret the results

Hope this isn't too much.   Sounds like a homework problem.

Ced


If you are really interested in voice recognition. Check out Hidden markovs model. It is used by several softwares to recognize your words


Does a matched filter give a gain in SNR?

Output SNR for matched filter can be analysed in time domain, which is as explained below,

Consider a filter with Impulse response  $ \enspace h(t)$  and input $\,s(t)+w(t)$, where $\,s(t)$ represents the required signal and $\,w(t)$ represents AWGN noise.

Output signal is given by  $$y(t) = s(t) * h(t) + w(t) * h(t) $$

Writing signal component of output in inverse fourier form   $$ y'(t) = \int_{-\infty}^{\infty}  S(f) H(f)  \ e^{j2pit}\ df $$ where $S(f)$ and $H(f)$ are spectrums of input signal $s(t)$ and Frequency response of filter.

Output SNR of filter which is sampled at time instant $\,t=T $ is output power with signal component by noise power which can be return as

$$SNR  =  \frac{\displaystyle\left |\int_{-\infty}^{\infty}  S(f) H(f)  \ e^{j2piT}\ df \right|^2}{ \displaystyle\left|\int_{-\infty}^{\infty}  \frac{N_o}{2} |H(f)|^2  df \right|^2}$$

Note that numerator term is evaluated at time t = T, which is the sampling instant for Matched filter.

To maximise this SNR , we need to cancel the phase term .
if we substitute $$H(f) = S^*(f)  e^{-j2pif} $$
which is nothing but $$h(t) = s^*(T-t) $$

and simplify the above SNR expression, we will get
$$ SNR = \frac{\displaystyle\int_{-\infty}^{\infty} \left|S(f)\right|^2 df}{N_o/2} $$
which is nothing but ration of energy of signal $s(t)$ to value of power spectral density of AWGN noise.

In this time domain analysis you can see that there is no filter effect , and it also tells us what should be the filter impulse response to maximise the output SNR.


This is about signal processing in the first radar

Indeed, in the very first time of radar, the technology necessary to measure delay between transmited and received echoes wasn't avaible for EM waves.
I think that these first type of CW radars weren't used for ranging, only for detection.
Using only one sensor which emit and receive the CW signal you can't perform ranging. In a bistatic configuration, i.e. the transmitter and the receiver arn't colocated, you may perform rangin by imposing a constraint on the target's kinematic (for instance, constant speed) and integrating measurments over time.
Using a network of sensor (so called multistatic configuration) you may perform "real time" ranging without having to integrate measurments nor imposing a kinematic model on the target's motion.

Direct answer : These early radar were used only for detection and not for ranging.


Frequency of Signal of a point tracked in a video?

You could use Fourier transform by using FFT then find the magnitude of each frquency or simply find the power spectrum by functions like pwelch.


Signal corrupted using peak (biquad) filter at low frequencies

Looks a like a typical stability problem with the poles too close to the unit circle. A third octave wide parametric EQ filter at 25Hz sampled at 44.1kHz has a pole radius of > 0.999 so it's getting VERY close to the unit circle. The lower the frequency, the large the pole radius.The pole section alone has a peak gain for more than 110 dB !!!

With the poles so close to being unstable, you may run into all sorts of problem: numerical noise, noise amplification, coefficient rounding or truncation errors, filter topology (do NOT use Direct Form II or transposed Form I), etc.


No, it will be Rician Distribution.

Rayleigh distribution is a special case of Rician Distribution when the normal random variables involved are of zero mean and equal variances.

I have assumed you meant :$$y=x_r + j \cdot x_i + n_r + j \cdot n_i$$ where $x_r, x_i$ are deterministic variables.

Basically, you now have $y_r = (n_r + x_r),\ \sim \mathcal N(x_r, \sigma^2) $ and $y_i = (n_i + x_i), \ \sim \mathcal N(x_i, \sigma^2)$

Now, you can write $y = y_r +j \cdot y_i = r \cdot e^{j\phi}$, where $r,\phi$ are random variables and, $y_r = r \cdot \cos(\phi)$ and $y_i = r \cdot \sin(\phi)$.

Using method of transformations, Jacobian matrix can be given as follows:
$$\begin{pmatrix} y_r \\ y_i \end{pmatrix} = \begin{pmatrix} r \cdot \cos(\phi) \\ r \cdot \sin(\phi)\end{pmatrix}$$
$$\mathbf J = \begin{pmatrix} \cos(\phi)&&-r \cdot \sin(\phi) \\ \sin(\phi)&& r \cdot \cos(\phi)\end{pmatrix}$$
$$|\mathbf J| = r$$
Hence, you see that Jacoboan does not change by changing the mean of the normal random variables from 0 to non-zero.

Joint PDF of derived RVs $r, \phi$ will be given by the following :
$$f_{R,\Phi}(r, \phi) = f_{Y_r, Y_i}(r\cos(\phi), r\sin(\phi)) \cdot |\mathbf J |$$
$$f_{R,\Phi}(r, \phi) = \frac{r}{2\pi \sigma^2} \cdot \exp\left[-\frac{(r\cos(\phi) - x_r)^2+(r\sin(\phi)-x_i)^2)}{2\sigma^2}\right]$$
$$\forall r\in [0,\infty) \\ \phi \in[-\pi, \pi]$$

Solve for Marginal PDF of $R$ and $\Phi$ to see that it is not Rayleigh anymore. But a more generalized Rician Distribution.


Check the stability of the system, $y[n]=u[n]$

(Long Comment) Sorry guys, I did not get why that system created so many conceptual problems to you. For instance, consider the unforced system $y(n)=0$. If I define $\pi:\{0\}\times\mathbb Z\to\{0\}$, $\pi(y,n)=0$ then $(\{0\},\mathbb Z,\pi)$ defines a dynamical system on the state space $\{0\}$. In fact $\pi$ is continuous, and it satisfies the identity axiom ($\pi(x,0)=x$ for all $x\in\{0\}$) and the group axiom ($\pi(\pi(x,n_1),n_2) = \pi(x,n_1+n_2)$ for all $(x,n_1),(x,n_2)\in\{0\}\times\mathbb Z$).

I agree that this is a degenerate case. I guess I can extend the state space to $\mathbb R$ by defining the system in a different way. For instance, let
$$
x(n+1)=0\qquad (1)
$$
be a ordinary difference equation and define $\pi$ as the flow of (1). Then $(\mathbb R,\mathbb N,\pi)$ defines a dynamical system on $\mathbb R$. Take the forced system
\begin{align*}
x(n+1) &= 0\cdot x(n) + 0\cdot u(n)\\
y(n) &= 0\cdot x(n) + 1\cdot u(n)\qquad\qquad (2)
\end{align*}
That makes sense to me. It is a LTI system with zero relative degree. To check BIBO stability let take a reference input $\bar u:\mathbb N\to\mathbb R$ and a perturbation $\delta u:\mathbb N\to\mathbb R$. Let $u(n)=\bar u(n)+\delta u(n)$.
Let $\bar y:\mathbb N\to\mathbb R$ be the solution corresponding to the input $\bar u$ and let $y:\mathbb N\to\mathbb R$ be the one corresponding to the input $u$. Then
\begin{align*}
\|y(n)-\bar y(n)\| = \|u(n)-\bar u(n)\|=\|\delta u(n)\|
\end{align*}
So trivially, for any $\epsilon>0$ there exists $\delta_\epsilon>0$ such that
$\|\delta u(n)\|\le \delta_\epsilon\implies \|y(n)-\bar y(n)\|\le\epsilon$. You can see this also as a continuity property of the solutions. In fact, you can define the space $Y$ of all the bounded output trajectories $y:\mathbb N\to\mathbb R$ and the space $U$ of all the bounded input ones. Let endow $Y$ with the norm $\|y\|:=\sup_{n\in\mathbb N}|y(n)|$ and $U$ with the norm $\|u\|:=\sup_{n\in\mathbb N}|u(n)|$. Then the system (2) definest an operator $T:U\to Y$, and asking for BIBO stability is the same as asking for the continuity of $T$ in the sup norms.


What is the difference between an equalizer and a channel estimator?

In communication systems, transmitted signals are distorted by the physical medium (the channel) charactheristics.

The channel estimator tries to identify the transmission channel characteristics, by computing (estimating) various of its model parameters.

Then, the channel equalizer tries to restore to distorted (received) signals back into their originals via some sort of inverse filtering (equalization) on the received signals, which was designed according to the estimated channel parameters.

Both operations are approximate; depending on the application, they can be very effective or not so.


A magnitude-only spectrogram is not invertible. The spectral phase required for reconstruction has been lost.

A complex FFT (a 1D column of spectrogram data before the absolute magnitude is taken and colored) might be (within numerical limitations) invertible by using an IFFT, but that won't recover the original signal if the STFT is using windowed FFTs.


Why does the frequency in the DFT have to be an integer?

There are a few different ways to interpret the math of the DFT. The one I find most suitable to explain this and other properties of the DFT assumes that the continuous Fourier domain was sampled, which causes the time domain to become periodic (just like the signal was sampled, causing the Fourier domain to become periodic).

If your signal contains a sine wave, it has a single frequency. But if a non-integer number of periods of this sine wave is represented by the samples in your discrete signal, then there will be a discontinuity when replicating these samples. This discontinuity contains all frequencies.

The "frequency" parameter used in the DFT, often represented as k, is in relation to the length of the signal. If your signal shows exactly a single period of the sine wave, this sine wave has a frequency k=1. If you see exactly two periods, it has a frequency k=2. If you see one and a half periods, then k would be 1.5, which does not exist.

Draw one and a half period of a sine wave. Now copy that signal multiple times end to end:

x = cos(linspace(0,3*pi,1000));
plot(repmat(x,1,7))
set(gca,'xticklabel',0:7)

You now have a signal where the lowest frequency corresponds to the repetition of the original signal, the period length is the length of the original signal. This is k=1. The next frequency is twice this lowest frequency (period is half). You will notice if you study the graph above that there is nothing that repeats with a frequency of 1.5 times the lowest frequency (period is 1.5 times the original signal). Such a frequency does not exist because of how the signal was constructed.


Mapping the Value of a Sample in a 2D DFT to Cycles/Pixel

Well the Nyquist frequency is 0.5 cycles / pixel and the range of frequencies can be considered to range from -0.5 to +0.5 cycles / pixel. So if the dimension in a given axis is N pixels then each pixel is 1 / N cycles per pixel.

Depending on your particular DFT/FFT implementation you may have zero frequency in the middle of the image (M/2, N/2) or in the corner (0, 0). For the former case the spatial frequency at u, v will be:

(M/2 - u) / M, (N/2 - v) / N


The way I prefer to look at it is to think of the QPSK signal as a complex envelope; that is, as a baseband signal. In this representation, you have a single basis, but the amplitudes are allowed to be complex.

So, let $p(t)$ be a root-raised cosine pulse such that $$\int_{-\infty}^\infty p(t-lT_p)p(t-mT_p) = \begin{cases}1\textrm{, if $l=m$}\\ 0\textrm{, if $l\neq m$,}\end{cases}$$ where $R_p=1/T_p$ is the pulse rate. Now you can think of the QPSK signal $s(t)$ as $$s(t)=\sum_k a_kp(t-kT_p)$$ where $a_k$ are complex symbols taken from a QPSK constellation. As you can see, you only need one set of real orthonormal pulses. The complex amplitudes $a_k$ can be recovered by matched filtering.

Another way to look at it is that you have two basis, one real and the other imaginary (the set $jp(t)$).


What's wrong with my Goertzel algorithm implementation?

Why do you think the two should be equal?

Check out the final line of equation (6) in the article you link to:

I believe your misunderstanding is to not take account of the $e^{j\omega_0n}$ term.


Librosa stft + istft - Understanding my output (which always seems too perfect) at varying window lengths

The information you think is missing, due to either poor visible frequency resolution or poor visible time resolution, is actually still contained in the sequence of FFTs, but not in a form that can be easily visualized, especially by just displaying FFT magnitudes or spectrograms.  That's because the information is "hidden" in the FFT phase results.

In the sequence of short FFT windows case, some of the frequency resolution information is represented in the difference in phase between adjacent or overlapped windows.  In the long FFT window case, much of the time locality information is represented in both the absolute phases and in the difference in phases within groups of nearby FFT result bins.  None of that is visible in FFT magnitude plots or spectrograms; and FFT phase and phase delta plots are hard to read.

But complete complex IFFTs do use the phase information contained in the real and imaginary components of the complex FFT results, and thus can reproduce both the time and the frequency information accurately and completely.

This phase information can be used in certain forms of FFT analysis to gain back some of the "missing" resolution.  For instance, the phase vocoder algorithm deduces or uses greater frequency resolution by looking at the phase deltas between adjacent or overlapped FFT windows.  And the time shift property of the FT shows how the phases between DFT result bins rotate relative to each other as an event (an AM or "beat" modulation, etc.) is rotated to different positions within a DFT window.


Effect of nonlinear transformation on digital signal

Your intuition is correct. When simulating a discrete system, the sampling frequency needs to be high enough that you meet the Nyquist sampling criterion at all times, or you'll have folding.

A method I've used to determine the sampling frequency is to assume the signals are continuous, and find the Fourier transform of the system output. Usually you have an infinite number of harmonics, but you can usually identify a threshold after which you don't care about the harmonics folding. What the threshold is depends on your application and how much folding can you tolerate. Usually, harmonics 30 dB below your signal should pose no problems.
You can't suppress the harmonics at the nonlinear transformer's output, but you can filter them afterwards. Note that this will not help with the sampling frequency, since the harmonics are present between the transformer and the filter.
You need to make sure the signal is sampled at a frequency the DAC accepts. You can achieve that by resampling the signal right before the DAC.


Programming the IDWT for image processing

here is one version of the Haar analysis/synthesis filter bank pair. Oddly, the synthesis filters, up to a scale factor, are mirrors of  the analysis filters.

Other options are:

$\pm 1/\sqrt{2}$ instead or $1/2$ and $1$ (orthogonal version);
reverse: $[1 \;  \pm 1]$ at analysis and $[\pm 1/2 \; 1/2]$ at the synthesis.


Inverse Z-transform with complex conjugate poles

Since the given $\mathcal{Z}$-transform

$$X(z)=\frac{2+3z^{-1}}{1-z^{-1}+0.81z^{-2}},\quad |z|>0.9\tag{1}$$

has real-valued coefficients, its poles and zeros must be either real-valued or they must occur in complex conjugate pairs. The poles of $(1)$ are a complex conjugate pair:

$$p_1=0.5 + 0.74833j=re^{j\phi}=p\\
p_2=0.5 - 0.74833j=re^{-j\phi}=p^*$$

with $r=0.9$ and $\phi=0.98177$.

You can proceed from the partial fraction expansion of $X(z)$, which you correctly obtained:

$$X(z)=\frac{A}{1-pz^{-1}}+\frac{A^*}{1-p^*z^{-1}}$$

With the given region of convergence (ROC) $|z|>0.9$, we know that we're looking for a causal sequence $x[n]$:

$$x[n]=\left(Ap^n+A^*(p^*)^n\right) u[n]\tag{3}$$

With $p=re^{j\phi}$ and $A=|A|e^{j\phi_A}$, $(3)$ can be rewritten as

$$x[n]=|A|r^n\left(e^{j(n\phi+\phi_A)}+e^{-j(n\phi+\phi_A)}\right)u[n]=
2|A|r^n\cos(n\phi+\phi_A)u[n]\tag{4}$$

which is of course real-valued.

With the given values of $p$ and $A$ the impulse response looks like this:


I believe that h(-t) means a "time-reversed" version of h(t). Your command: 'y = conv(r,-h);' computes the convolution of 'r' and negative 'h', and you don't want that. I think you want:

y = conv(r,conj(fliplr(h)));


feature detectors and descriptors comparison

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


Impact of the WSOLA analysis increment on the time stretching factor

Yes, there can be jitter between the requested alpha and the one provided by WSOLA; but that jitter will average out over the long run, and the peak jitter is bounded by on the order of twice the lowest frequency (half the longest period) accepted by the autocorrelation algorithm.


Question on the update step of Bayesian filter

I would prefer to see the definitions, but based on your prior work (which I haven't checked thoroughly), we can already reach the final stage. You left off at $$ p(x_k | y_{1:k}) = \frac{p(y_k|x_k) p(x_k|y_{1:k-1})}{p(y_k|y_{1:k-1})} $$ which is a probability distribution of $x_k$ based on the left hand side. The term in the denominator is independent of $x_k$ (if I have guessed your model correctly). Therefore, it is an unknown positive constant. That is, there exists some constant $Z_k$ such that $$ p(x_k | y_{1:k}) = \frac{1}{Z_k} p(y_k|x_k) p(x_k|y_{1:k-1})$$ where $Z_k = p(y_k|y_{1:k-1})$.

Now let's integrate both sides.
$$ \int p(x_k | y_{1:k}) d x_k = \int \frac{1}{Z_k} p(y_k|x_k) p(x_k|y_{1:k-1}) d x_k $$
The left hand side is a probability distribution with respect to $x_k$, therefore it sums to $1$. Because $Z_k$ is independent of $x_k$, we can pull it out of the integral. Multiply both sides by $Z_k$ to obtain your answer,
$$Z_k = \int p(y_k|x_k) p(x_k | y_{1:k-1}) d x_k.$$

Is $Z_k$ a constant with respect to $x_k$?

We know that $Z_k = p(y_k|y_{1:k-1})$ in other words $Z_k$ is the probability of $y_k$ given all previous $y_{1:k-1}$ after integrating out $x_k$. Note that this is the same thing as $$Z_k = \int p(y_k|y_{1:k-1},x_k) p(x_k) d x_k$$
which (after integrating) removes all dependency on $x_k$. Therefore, as noted in the original response, $Z_k$ is constant with respect to $x_k$ and can be factored out of the integral.

This practice of marginalizing out a variable based on its prior is a common theme in Bayesian systems.


The consumption time and transmission time is identical: One second of data is still one second of data regardless of sampling time.

The greatest common divisor between the two rates is 300, thus to resample this exactly from 44.1KHz to 48KHz you would need to use the ratio $160/147$ (and the inverse for the other direction):

$147$ is factored into $3, 7^2$

$160$ is factored into $2^5, 5$

The following demonstrates one approach to resample from 44.1KHz to 48KHz, where care has been taken to not reduce the sampling rate below 44.1KHz (if that matters for fidelity concerns) and the multiple stages simplifies the filtering needed:

Interp by 4, decimate by 3, interp by 8, decimate by 7, interp by 5, decimate by 7.

This would be implemented with the following structure where the interpolator blocks signify insert of $I$ samples between each sample and the decimation blocks signify selecting every $D$th sample and throwing away the rest. The intermediate blocks can run at any arbitrary higher sampling rate to keep up with the throughput and the input/output blocks are rate matched (consuming samples at 44.1KSps and providing output samples at the 48 KSps rate).  To do this as shown with 20 MHz audio bandwidth and 80 dB resampling image rejection, I estimate that 171 taps would be needed for FIR1, 95 taps for FIR2 and 25 taps for FIR3 (as linear phase filters so one multiplier for every 2 taps).  For real-time application, the expected delay through the resampler would be 7.9 ms.

The filters could certainly be designed with windowed Sinc functions (this is known as the windowing approach to FIR filter design which is sub-optimal - see our further discussion here FIR Filter Design: Window vs Parks McClellan and Least Squares). The least squares algorithm (firls in MATLAB/Octave and Python) provides an optimal solution for resampling applications, resulting in higher image rejection for a given number of taps. Further, in many resamplers (not this one due to the close ratio), the images to be rejected can be isolated to distinct frequency bands; resulting in the use of multiband filters which the least squares algorithms support and further maximize rejection where it is needed most.

Interpolation and decimation resampling can also be accomplished by mapping the same filter coefficients as designed for the resampler shown above into polyphase structures as depicted in the diagram below. This would be identical in performance to the resampler above but can be done with an internal sampling rate as low as 67.2 KSps and significantly fewer overall computations. The decimation is done by selecting the appropriate filter output associated with the computation cycle for each decimator output rate (for the first two stages this just ends up being that the commutators move back one sample after every output update, and the last stage moves forward one sample after every other update). This can be a very efficient approach since only one of the filters in each stage actually needs to be computed for each output (note that each filter within a group would contain the exact same data but the multiply and sum only needs to be done on one of them each time). Since only one filter in each stage actually needs to be computed on any given cycle, the implementation can be done with just three FIR filters (44-tap, 12-tap and 5 tap) with the coefficients updated from a ROM table for each decimator output computation cycle (this high efficiency approach would require a tightly synchronized state machine, while instead computing all internal filters would allow for a lot of slop with the internal timing above the minimum limits). If 80 dB of image rejection or 20 MHz of audio bandwidth is not necessary, then all filter lengths can be reduced accordingly.

Resampling with polyphase filters are detailed further in this post:

How to implement Polyphase filter?

As Robert points out in the comments, the above would work continuously if the input and output were at the exact frequencies given; or if the input frequency was slightly less as the processing can ensure the next sample is ready prior to being clocked by the output. The problem will be if the input frequency error is slightly higher (or output frequency slightly lower) that lost samples will occur. Some buffering can be provided to sustain modest frequency variation but will invariably overflow. The only robust solution for a real time application without operating time limitation is to ensure the input and output clocks are synchronized through some mechanism.  Given the OP is mentioning sending the data over the network, I believe buffering would then be required based on the predicted worst case clock inaccuracies between the two locations and longest time duration of an audio transmission.  If the input and output clocks were co-located then they could be PLL locked to each other.


Firstly the terminology: The thing that you are trying to compute is the Discrete Fourier Series (DFS). This the only flavor for Fourier Transforms that is discrete in both time and frequency and can thus be represented numerically inside a computer. The Fast Fourier Transform is a specific algorithm to computer the DFS (see for example http://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm ). The term FFT is often used to refer to the DFS but that's actually somewhat wrong and often confusing.

You can always calculate the DFS directly using it's definition but this will require $N^{2}$ complex multiplies. The FFT speeds this up by breaking splitting the vector N into smaller sub vectors and then using the symmetry properties of the transform coefficient (often referred to as "twiddle factors"). The vector needs to be split in equal pieces though and that works better the more prime factors N has. A power of two is the best since it has the most prime factors and the factors themselves are the smallest possible. Worst case are N's that are prime themselves.

Let's assume you want to calculate the DFS of sequence of 1021 points. Since it's a prime number you need to use the direct formula which requires 1021*1021 multiples, roughly a million. If you zero pad to 1024, you can use the most efficient FFT version and the number of complex multiplies is 2*log2(N)*N is only about 20000. Hence that's a lot faster.

FFT is efficient if your DFS length can be broken down into lots of small prime numbers. Large prime numbers are slower. In practice you can typically find a spot made our of 2s and 3s that's fairly close to where you want to be so there is no need to always go to a power of 2.


blob detection using Laplacian of Gaussian filter

Most likely your images look different from the ones in the lectures because of scaling. Note that the result of the convolution with a Laplacian filter will have positive and negative values. What the resulting image looks like depends on the data type of the array, and on the range to which the values are scaled.  For example, if you store your filtered image in uint8, without doing any scaling yourself, the negative values will be truncated to 0.

Typically, when an image has negative values, it is scaled so that 0 corresponds to a neutral gray, the minimum value to black, and the maximum to white. It looks like this is what is happening in your code.

On the other hand, in this particular case you are probably interested in both maxima and minima of the Laplacian (dark on bright and bright on dark blobs), so you may as well just take the absolute value of the filtered image.  Then you will not have any negative values, which makes the scaling easier.

As far as your question about the edges, this is exactly what should be happening. The value of the Laplacian at an edge should be close to zero. This property is used for very precise edge localization by finding the zero-crossings of the Laplacian of Gaussian.


problem finding 16 point DFT using two 8 point FFT (Divide and combine algorithm) MATLAB

I will only give insight how to work on above problem but I won't give you the matlab code it can solved using two methods.
Method_1:- using DIT FFT algorithm hear inputs given in bit reversal order and output is multiplied with twiddle factor$w_n^k,where  k=0,1,2..,7$

So you use the same for your problem.
Method_2:- using DIF FFT algorithm hear the inputs are given in normal order at the output the butterfly's used are 2X2 & output is taken in bit reversal order like below


I'm not sure what you're getting at with "reflection" or "striping," but you want to duplicate your data points such that the frequency-domain signal (the input to the IFFT) is conjugate symmetric. That means, for $N$ even:

$$
X[k] = X^*[N-k], k = 1, 2, \ldots, \frac{N}{2}-1
$$

Note that $k$ is zero-based; the first element in the vector is $X[0]$. The above indicates that you do not duplicate the zeroth and $\frac{N}{2}$-th elements in the input vector. This relationship is based on the periodicity of the discrete Fourier transform; the first half of the frequency-domain vector corresponds to "positive" frequencies covering the angular frequency range from $0$ to $\frac{(N-1)\pi}{N}$. The second half covers from $\pi$ to $\frac{(N-1)2\pi}{N}$, which, taking advantage of the $2\pi$ periodicity in the frequency domain, is equivalent to $-\pi$ to $\frac{-2\pi}{N}$. In order for the frequency-domain vector to have conjugate symmetry, then it can be easily shown that the above relationship must be true.

Speaking to your original goal of converting a spectrogram into an audio signal that sounds meaningful, there's no guarantee that what you get out will be in any way pleasing. As was pointed out on Stack Overflow, the spectrogram has no phase information; the method suggested there assumes the phase of each frequency-domain bin is zero (so the conjugate operation shown above is superfluous). While the human ear isn't terribly sensitive to phase distortion in audio, you might have to do something else depending upon exactly what you're hoping to achieve.


Read the book 4G Lte Advanced by Stefan Parkvall and Eric Dahlman. According to the book,
FFT size of 2048, with a corresponding sampling rate of 30.72 MHz, is suitable for the wider LTE carrier bandwidths, such as bandwidths of the order of 15 MHz and above. However, for smaller carrier bandwidths, a smaller FFT size and a correspondingly lower sampling rate can very well be used.(Lowering the sample rate simplifies the implementation in dual-mode HSPA/LTE terminals)

Also,nothing forbids the implementation of an OFDM modulator as a set of parallel modulators,
Also, nothing prevents the use of a larger IFFT size, for example a size-2048 IFFT size, even in the case of a smaller number of OFDM subcarriers.

Tcp= 160·Ts and Tcpe=512.Ts. Change in CP induces change in the total number of REs(6X12, instead of 7X12). Note that #sub carriers remain 12 in both cases and FFT size is dependent on BW allocated which is expressed in #sub carriers.


Manually creating a lowpass filter from the frequency domain: what phase?

Typically if you're designing to an amplitude response, you set the phase at zero.  In processing time-domain signals, this is to keep the group delay of the filter constant.  In processing images, this will keep the filter effects centered, rather than displacing different spectral components of the image with respect to one another.
At least in time-domain processing there can be good reasons not to do this.  I can't think of any reason you would want a non-constant delay filter in image processing.
Keeping the phase at zero makes a constant-delay filter.  If you're doing your filtering after the fact, then you can just delay everything else by the same amount so that it appears that you're using a zero-delay noncausal filter.
However, if you're responding to the filter output in the real world (usually for some sort of closed-loop control, but possibly for fault detection), then that delay can be a performance-killer.  In addition, there's a feeling in the audio community that unwisely-designed constant-delay filters can cause aesthetic problems (Google on "pre-ring").  You solve these problems (and cause others) by using a minimum-phase filter.  These are almost always done by using an IIR filter rather than an FIR.


How to convert sampled real audio signal to it's IQ components for IQ modulation?

I’ll answer the question I think you’re asking here:

If you have a real valued signal that came right out of an A/D and you’d like to recover signal phase, you can use the Hilbert transform to generate a version of the signal that is in quadrature with your A/D output. In this case, your original real valued signal would be the I component (in-phase) and the output of the Hilbert transform multiplied by j, the imaginary number, would be your Q component (quadrature).

Now if you have some modulation you want to deal with, you would do the IQ recovery first like I just described; once you have that IQ signal you’re pretty much completely characterized the signal, so you can do just about whatever you’d like with it while preserving phase.


What do we mean by Total Factor Space in Audio Processing?

I read the paper, let me try to summarize what I understand.

1) For speaker identification tasks, systems need to be robust to  "session variability". This means that from one spoken utterance to another, the way a speaker pronounces the same passphrase & the environmental/background noise impacts adds a bit of variance to the captured audio signal. I.e. it is not exact.

2) Historically, speaker-verification techniques based on GMM-UBM->MAP are not robust to channel distortion (i.e distortion during the recording) or phonetic variability (i.e. differences in pronunciation) or other "nuisance factors"

3) A single space named Total Factor Space is constructed to model speaker and channel variability jointly. Joint factor analysis (JFA) is one technique which decomposes a speaker's "mean supervector"  (which is a linear combination of the set of eigenvoices) into speaker component S and channel component C

So instead of doing separate modeling of speaker and channel spaces, all factors (i.e. both speaker and channel) are combined into one space. The JFA model from the slideshare below looks like:

Have a deeper look at this slideshare. Good luck!


FFT question as it relates to Neural Networks / Supervised Learning Models

I stumbled across this.
http://ataspinar.com/2018/12/21/a-guide-for-using-the-wavelet-transform-in-machine-learning/

hope it helps
john


What are the problems with designing an FIR filter using FFT?

One reason you see people designing FIR filters, rather than taking a direct approach (like both 1 and 2) is that the direct approach usually fails to take into account the periodicity in the frequency domain, and the fact that convolution implemented using an FFT is circular convolution.

What does this mean?

Suppose you have a signal $x = [ 1, 2, 3, 4]$ and a filter impulse response (convolution kernel; you are correct they are the same) $h = [ 1, 1 ]$.

The convolution $y = x *h$ is $[1, 3, 5, 7, 4]$, a 5-length vector.  If you use the FFT (of the wrong length, 4) then the answer you get is $[3, 5, 7, 5]$.  The reason for the difference is that the result of linear convolution of these two is length 5, but the result of circular convolution is whatever the FFT length was.

If the FFT length is greater than or equal to the length of the result of linear convolution, then the two are the same.  Otherwise, the two are not the same (unless the data somehow conspires to make it so, e.g. if one signal was zero).


For implementing the frequency shift of 7.5 kHz you will have to apply oversampling at some point in your system. There are two principal methods to achieve this: in frequency domain or in time domain.

Frequency domain

You could simply apply an additional zero padding in frequency domain by using a 2048 IFFT and by setting every second subcarrier to zero. You then shift all subcarriers by one frequency bin to achieve the frequency shift. This will maintain the frequency spacing of 15 kHz but the bandwidth of each subcarrier will be divided by two compared to the 1024 IFFT.

Time domain

You apply two-times oversampling to the time domain signal at the output of the IFFT. This can be done by padding a zero sample after each sample followed by a lowpass filter or some other interpolation method. The frequency shift can then be achieved by multiplying the resulting signal by $\exp(j\Omega_1 n)$, where $\Omega_1$ is the 7.5 kHz normalized to the new sampling frequency (double the original sampling frequency) and $n$ is the discrete time.

Which method to choose depends on how close your simulation should be to the hardware implemenation. If the system you model has some intermediate frequency oscillator to do the frequency shift, the time domain method can model this. But if your digital-to-analog converter is fast enough and if you can implement a 2048-IFFT at the required speed, this is probably the way to go and you should therefore model it accordingly.


How to find the value of regularization parameter/s when having one or multiple priors

I have just found this paper that calculates the $\lambda_i$ parameters on an L-surface: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.135.5969&rep=rep1&type=pdf


How to evaluate CCDE at a given frequency?

You need to follow the following steps (reponses at each steps are in a spoiler block, try to find them by yourself before checking). I will use the usual notation $z=e^{j w}$ for simplicity.

Find the relation between the transfer function $H(z)$ and the discrete-time Fourier Transform of the input $X(z)$ and the output $Y(z)$. This should be a basic results of your course

   $H(z) = \dfrac{Y(z)}{X(z)}$

Take the (discrete-time) Fourier Transform of your CCDE. For that, you will need to use two basic properties of the TF: linearity and time-shifting

   $Y(z) - \frac{1}{2} z^{-1} Y(z) = 2 X(z) - 5 z^{-1} X(z) - z^{-2} X(z)$

Compute $H(z)$ from the previous equation

   $H(z) = \dfrac{2 - 5z^{-1} - z^{-2}}{1-\frac{1}{2} z^{-1}}$

Finally take the value for $z = e^{j \pi} = -1$

   $H(e^{j \pi}) = \dfrac{2 + 5 - 1}{1+\frac{1}{2}} = \dfrac{6}{\frac{3}{2}} = 4$

This last result mean that your filter, at Nyquist frequency, has

   no phase delay and a gain of 4.


Difference Between Equiripple & Least Squares Design for FIR Digital Filters

While I completely agree with Jason R's answer, I would like to add a few things that I consider important. First of all, it is a misunderstanding to believe that for least squares designs the transition band width is smaller than for equi-ripple designs. The width of the transition band depends on many design parameters but it is independent of the optimality criterion. With both criteria you can define arbitrary transition bands ("don't care regions") where you simply do not specify a desired response. This has to be done explicitly for the Parks-McClellan algorithm, but it can (and should) also be done for least squares designs. The most basic least squares design which is to simply truncate the Fourier series of a (often discontinuous) desired frequency response, is definitely no benchmark for comparing least squares designs with other optimality criteria.

The second important thing to notice is that the least-squares design and the equi-ripple design are two extreme points on a trade-off curve between maximum error and error energy. While the least squares design minimizes the error energy, its maximum error is relatively large, and the opposite holds for the equi-ripple design. So it is desirable - and possible - to mix the two criteria. One way of doing this is to use a constrained least squares criterion which minimizes the error energy while constraining the maximum error to some desired limit. In many cases it is desirable to have an equi-ripple design in the passband of a filter and to minimize the (weighted) least squares error in the stopbands, because then the maximum distortion of the signal is minimized, while the noise power in the stopbands is also minimized. The following example shows such a design. It is a linear-phase FIR bandpass filter of order 100, with a constraint on the maximum error in the passband (in this case resulting in an equi-ripple design in the passband), and with minimal error energy in the stopbands. The top figure shows the magnitude of the frequency response in dB, and the bottom figure shows the magnitude of the frequency response in the passband. Clearly, we have an equi-ripple approximation in the passband, and a least squares approximation in the stopbands.


Hilbert Transform of a sine Function with Quadratic Argument $\sin(At^2 + Bt + \pi/4)$

Directly calculating this seems difficult.

My argumentation is the following. For a signal $s(t)$, the so-called analytic signal $s_{\mathrm{an}}(t)$ can be obtained by
\begin{equation}
s_{\mathrm{an}}(t) = s(t)+\mathrm{j} \mathcal{H}\{s(t)\}\,\,, \text{where}\,S_{\mathrm{an}}(\omega) = 0\,\forall\, \omega < 0
\end{equation}
The analytic signal essentially corresponds to the spectral content of $s(t)$ in the positive frequencies only.

For your first example of a simple sinusoid, you can also come to the result for the Hilbert transform if you consider the analytic signal. The real sinusoid consists of the frequency components $\pm B$. The analytic signal should then be the component at $B$, which is then obviously a complex exponential, yielding the result for the Hilbert transform.

Now for your chirp signal $x(t)$, the situation is somewhat more complicated. If we think about a virtual "instantaneous frequency course" of the signal, it is
\begin{equation}
\omega_{x}(t) = \pm (2At+B)\,.
\end{equation}
Now this is somewhat strange, corresponding to two linearly changing components of opposite slope, crossing the zero-frequency point at $\omega_{x}(t=-\frac{B}{2A})=0$.

Now the analytic signal would have to represent the part of this frequency course above the zero-line of the $\omega-t$ plane (I might add some plots later). This means it would have to first have a negative slope, go down to frequency zero, and then abruptly change to a positive slope!

This means that the analytic signal would have to look something like
\begin{equation}
x_{\mathrm{an}}(t) = c_1 \exp{(-\mathrm{j} (At^2 + Bt + \frac{\pi}{4}))}\,\forall\,t < -\frac{B}{2A}
\end{equation}
and
\begin{equation}
x_{\mathrm{an}}(t) = c_2 \exp{(\mathrm{j} (At^2 + Bt + \frac{\pi}{4}))}\,\forall\,t \geq -\frac{B}{2A}\,\,,
\end{equation}
where the $c$ are some constant with $\vert c \vert = 1$.

Now we can determine the Hilbert transform of $x(t)$ by observing and checking in the equation for the analytic signal. This yields
\begin{equation}
\mathcal{H}\{x(t)\} = \cos{(At^2 + Bt + \frac{\pi}{4})}\,\forall\,t < -\frac{B}{2A}\,,\,\text{with}\,c_1 = -\mathrm{j}\,,
\end{equation}
and
\begin{equation}
\mathcal{H}\{x(t)\} = -\cos{(At^2 + Bt + \frac{\pi}{4})}\,\forall\,t \geq -\frac{B}{2A}\,,\,\text{with}\,c_2 = \mathrm{j}\,.
\end{equation}

One can probably also write these as one equation with the absolute value function. In any case, the point is the Hilbert transform seems to contain a discontinuity, which is what makes this particularly confusing to calculate I suspect.

I know it is somewhat "handwaivy", but I think the general idea/result is correct, so hope this helps!


You can always use normalised cross correlation with a structuring element like a disk or ellipsis to account for different viewing angles and marker distances. This is basically template matching, where a target that "looks like" the template is searched for in a larger image.

If more computational power is available, you could go down the route of detecting circles via the Hough Transform.

In this case, the Hough transform is used to identify some strong features that appear as points in some space and then, based on assumptions about their locations, the position of the circles is inferred. Similar to those techniques, are others that identify generic features off of an image (e.g. corners) and then try to cross-reference the position of the identified features with the relative positions identified off of a template.

A compromise might be to use similar methods as used to detect shotCodes. Now, I did try to locate that paper by the original authors of shotCodes where the detection process was outlined in detail but I cannot seem to find it. Instead, here is one that talks about detection (and the (more recent) associated thesis). This was quick enough to be implemented on early smartphones. Perhaps it offers a compromise in terms of ease of detection and position accuracy off of images.

Hope this helps.


Frequency response of a filter - fft vs. freqz

If you have the filter coefficients, freqz is the way to go - the various optional input arguments and the expected outputs are documented in MATLAB which you can see using help freqz. Snippets given below:

" freqz Digital filter frequency response.
    [H,W] = freqz(B,A,N) returns the N-point complex frequency response
    vector H and the N-point frequency vector W in radians/sample of
    the filter:

[H,F] = freqz(B,A,N,Fs) and [H,F] = freqz(B,A,N,'whole',Fs) return
    frequency vector F (in Hz), where Fs is the sampling frequency (in Hz).  "

If you have filter's impulse response - which will be different from the filter coefficients - you can use it to find the frequency response of the filter using the FFT.


Need some references on sound quality evaluation and syncing

From your comment I take it you are looking for speech quality measures. As far as I know, there exist a couple of objective measures, along with subjective procedures, since quality of speech is kind of a vague concept which may not be the same to everyone.

Since your application seems to require automatic evaluation, there are a number of measures you might want to check out:

ITU-T Rec. P 835 for assessment of listening speech quality (no available reference signal)
The PESQ for listening speech quality, mentioned by Florian already (when the reference signal is available)
The PEAQ for assessment of audio quality (when the reference signal is available)
The E-model für conversational speech quality (network parameters available)

All of these are listed in this book, along with a range of subjective measures.
It seems to me that for your application, the ITU-standard is most reasonable since no clean reference signal is required.

Edit: What also may be of interest of you is the frequency-weighted segmental SNR (fwsegSNR), which incorporates the articulation index (AI). This measure apparently correlates well with subjective measures.

An overview of objective speech quality measures is also given in this Loizou paper available online.


Why is my filter unstable and self-oscillating in this case?

The question is not very specific so this will also be just a general answer.

The stability of a composite linear-time-invariant (LTI) system composed of smaller LTI systems cannot be deduced from the stability of the component systems if the composition introduces feedback. Instead, you should test the stability of the full composite system. The system is stable if the poles of the transfer function lie inside the $z$-plane unit circle. So write out the transfer function of the component systems and also the transfer function of the composite system in terms of the component transfer functions. Find the poles of the transfer function of the composite system. Write an inequality that is the stability condition, and you may be able to solve from that an inequality for the acceptable values of a system parameter of interest, in terms other parameters.


Frequency Axis of Discrete Fourier Transform (DFT) with Odd Number of Data Points

This can be answered simply by considering the definition of the $N$-point DFT:

$$X_N[n] = \sum_{k=0}^{N-1} x[k]e^{-j2\pi \frac {n}{N}k }$$

where it's easy to see that the DFT just compares your $N$-point input signal $x[k]$ to a sinusoid of frequency $\frac nN$.

Thus, the lowest frequency is always $0$, and the resolution is always $\frac{f_\text{sample}}{N}$, no matter whether $N$ is even or odd. That also means that the lowest non-DC frequency is also $\frac{f_\text{sample}}{N}$.


The answer is yes but one has to specify $B_n$ properly to avoide possible confusions. In case if one uses a pulse compression, the bandwidth through which the receiver collects the noise will normally be $B_n = \beta_c$. Then, the "new" signal-to-noise ratio should be written as:

$SNR = \dfrac{P_TG_TG_R\lambda^2\sigma{P_g}}{(4\pi)^3R^4(kT_{sys}\beta_c)} = \dfrac{P_TG_TG_R\lambda^2\sigma{\tau_p\beta_c}}{(4\pi)^3R^4(kT_{sys}\beta_c)} = \dfrac{P_TG_TG_R\lambda^2\sigma{\tau_p}}{(4\pi)^3R^4(kT_{sys})}$

As you see, the signal-to-noise ratio remained the same as for the pulse without modulation (for which $B_n = 1/\tau_p$).

So, the pulse compression helps to increase the range resolution but the signal-to-noise ratio remains the same as for the unmodulated pulse. This holds true if you do not change $\tau_p$ and the peak power $P_T$ when you modulate your pulse.

The confusion (and as you mentioned "conflicting opinions") comes from the fact that it is quite often forgotten that if you increase the bandwidth of your pulse then you also increase the bandwidth through which the receiver collects the noise since radar's receiver is matched to the pulse and therefore has a bandwidth $\beta_c$.

Your last equation is correct. But since you put $F_n$ - the noise figure- then it should be specified that this is the signal-to-noise raito at the output of the receiver that has a noise figure $F_n$.


Why is the first value in digital first-order IIR filter impulse response not the largest?

A causal first-order IIR filter is characterized by the following difference equation:

$$y[n]=b_0x[n]+b_1x[n-1]-a_1y[n-1]\tag{1}$$

with $x[n]$ the input signal, and $y[n]$ the output signal. The impulse response of that system can be computed via the $\mathcal{Z}$-transform or otherwise, and it turns out to be

$$h[n]=b_0\delta[n]+(-a_1)^{n-1}(b_1-b_0a_1)u[n-1]\tag{2}$$

where $\delta[n]$ is the unit impulse, and $u[n]$ is the unit step sequence. Assuming that the filter is stable, we can conclude that $|a_1|<1$ is satisfied. Consequently, $|h[0]|>|h[n]|$ for $n>0$ is satisfied if and only if

$$|b_0|>|b_1-b_0a_1|\tag{3}$$

holds.

With $b_0=b_1$ and $a_1<0$ as in your example, condition $(3)$ can never be satisfied, and, consequently, the maximum magnitude of the impulse response is not assumed at index $n=0$.

EDIT:

I add more details referring to your filter design problem. A Butterworth low pass filter's frequency response satisfies $H(1)=1$ (DC) and $H(-1)=0$ (Nyquist). With

$$H(z)=\frac{b_0+b_1z^{-1}}{1+a_1z^{-1}}\tag{4}$$

we have the following requirements on the coefficients:

$$H(1)=\frac{b_0+b_1}{1+a_1}=1\\
H(-1)=\frac{b_0-b_1}{1-a_1}=0\tag{5}$$

which gives

$$b_0=b_1\quad\textrm{ and }\quad a_1=2b_0-1\tag{6}$$

From the stability condition $|a_1|<1$ we get

$$|2b_0-1|<1\Longrightarrow 0<b_0<1\tag{7}$$

Combining the constraints $(6)$ and $(7)$ with the inequality $(3)$ yields

$$1 >b_0>\frac12\tag{8}$$

If $(8)$ is satisfied, $h[0]$ will have the largest magnitude, otherwise it won't. Now it turns out that the exact value of $b_0$ depends on the desired cut-off frequency of the low pass filter. For cut-off frequencies less than $\frac12$ (i.e., less than half the Nyquist frequency), you get $b_0<\frac12$, and hence $|h[0]|$ is smaller than the maximum magnitude of the impulse response. Only for cut-off frequencies greater than half the Nyquist frequency do you get $|h[0]|>|h[n]|$, $n>0$.

So to directly answer the question in your title: the first impulse response value is not the largest (in magnitude) because the filter coefficients do not satisfy inequality $(3)$, and in your example of a (first-order) Butterworth low pass filter this happens because the cut-off frequency is smaller than half the Nyquist frequency.


Why is this MATLAB butterworth filter giving me wrong results

Your filter is passing only between 1 and 63 Hz, assuming that you meant to say that the sampling rate is 128Hz.

Since you are removing everything below 1Hz, the plot you have looks reasonable.

The range of values in your original signal is relatively small (4395 to 4420,) and the range of values in the output plot is also relatively small (cannot tell exactly due to scaling, but it looks small.)

The high values (+4000, -2000) in your output plot are the step response of your filter.  Your data jumps from zero (before the first data point comes in) to around 4400 at the beginning.  This step causes the filter to overreact - the short peaks up to 4000 and down to -2000 are the result.  This damps out right away, and the rest looks pretty much like what would be expected.

The offset (about 4400) looks like a DC component to the filter.  Since DC is lower than 1Hz, it gets removed as well.


I completely agree with Juancho's answer. I would just like to add a few things. I think the main problem is the misunderstanding which becomes evident in the last sentence of the question: "... can the Dirac function be sampled?" The Dirac impulse is NOT an ordinary function having definite values for every $t$, but it is a distribution (even though it is often called 'Dirac function'). One should therefore not try to 'evaluate' (or sample!) it. What is important about the Dirac impulse are its integral properties:

$$\int_{-\infty}^{\infty}\delta(t)dt = 1$$
and
$$\int_{-\infty}^{\infty}\delta(t-t_0)f(t)dt = f(t_0)$$

As Juancho already pointed out, the square of a Dirac impulse $\delta^2(t)$ is not defined. So if you were to sample a Dirac impulse, you would get the undefined result

$$\sum_n\delta(t-nT)\delta(t)=\delta^2(t)$$

Dirac impulses are a convenient tool for analyzing linear time-invariant systems but they should be treated with care because common types of  processing performed on ordinary signals (such as sampling) may lead to undefined and meaningless results when applied to Dirac impulses.


What is the meaning of the DTFT of the unit impulse sequence?

Your calculation is correct.

Your signal has its energy equally distributed in all frequencies.

But note that in discrete-time, "all frequencies" means from $\omega = -\pi$ to $\omega = +\pi$. Frequencies beyond that range exist but result in the same signal as one with a frequency within the range.

For example, $\cos(0\times n) = \cos(2\pi \times n)$. Similarly $\cos(0.1 \pi \times n) = \cos(2.1 \pi \times n)$, and so on.

So the energy of the signal is calculated integrating in a range of $2\pi$.

The formal result is known as Parseval's theorem, which states that energy measured in time domain is equal to energy measured in frequency domain:

$\sum_n |x[n]|^2 = \frac{1}{2\pi} \int_{-\pi}^{\pi} |X(e^{j\omega})|^2 d\omega$

Note: this is a particular version of Parseval. The more general version states that the projection (think cross product) of one vector onto another is the same calculated in time domain or frequency domain, which is reasonable since Fourier can be seen as a change of base, and a change of base should not alter projections.


Well, first of all you should note that you are taking batches of 512 samples. It basically means that the first FFT will perform a transform of the first 512 samples, slightly more than a quarter of a period, assuming that signal is periodic. (The windowing tries to compensate for this to some extent.) However, the information from that signal is not really what brings you any information about the spectrum.

While I haven't read the code, I suggest that you try one of these things (maybe both):

Increase the signal frequency
Increase the window length

Ideally, there should be at least one full period of the sinusoidal in the window if you want to interpret the result in a straightforward way.


Beginner Question: Extracting Features for Image Classification

You are just using the histogram as your features, yet your image contains very slight differences. This is not discriminative enough for SVM, you should focus on better features. In this case, you might be better of with subtraction, change detection and etc.


We have
$$G(z)=\dfrac{z^{-1}}{1-e^{-3}z^{-1}}$$
and
$$H(z)=\dfrac{1}{1-z^{-1}}$$
and we want to calculate
$$\dfrac{C(z)}{R(z)}=\dfrac{G(z)}{1+H(z)G(z)}$$

So let's do that:
$$
\require{cancel}
\dfrac{\dfrac{z^{-1}}{1-e^{-3}z^{-1}}}{1+\dfrac{1}{1-z^{-1}}\dfrac{z^{-1}}{1-e^{-3}z^{-1}}} = \frac{z^{-1}(1-z^{-1})}{(1-z^{-1})(1-e^{-3}z^{-1}) + z^{-1}}\\
= \frac{z^{-1}(1-z^{-1})}{1 - z^{-1} - e^{-3}z^{-1} + e^{-3}z^{-2} + z^{-1}}
= \frac{z^{-1}(1-z^{-1})}{1 \cancel{- z^{-1}} - e^{-3}z^{-1} + e^{-3}z^{-2} \cancel{+ z^{-1}}}\\
= \frac{z^{-1}(1-z^{-1})}{1 - e^{-3}z^{-1} + e^{-3}z^{-2} }
$$
so it looks like either the solution you have is wrong, or that $G(z)$ is not right.

So, let's see what could be wrong. Either
$$
z \not = e^{sT}\tag{1}
$$
or
$$
s \not = 1 - z^{-1}\tag{2}
$$
or
$$
G(z) \not = \dfrac{1}{1-e^{-3}z^{-1}} \cdot {\cal T}_{s \rightarrow z}(e^{-s})\tag{3}
$$
where ${\cal T}_{s \rightarrow z}$ is the $s$ to $z$ conversion.

Comparing (1) and (2), it looks like (1) is assuming the starred transform and (2) is assuming the impulse invariance method --- which can't both be right.


Help in understanding the formula of Signal-to-Noise-Ratio (SNR) - Part 1

$\mathsf{SNR}$ (signal-to-noise ratio) is a generic term whose value can be defined in different ways by different people, and as long as one states clearly what is meant by $\mathsf{SNR}$ in a particular document, there is no confusion. Thus, there is no "arriving" at the formula
$\mathsf{SNR} = \frac{E[y^2]}{\sigma_w^2}$ at all: it is the definition
of the term $\mathsf{SNR}$ as it used in that particular document.  There is, of course, a general feeling among communications systems designers and analysts that the error probability $P_e$ should be a decreasing function of $\mathsf{SNR}$, but exactly which decreasing function it is depends on the definition of $\mathsf{SNR}$ in use. For example, someone could define $\mathsf{SNR}$ as $\frac{E[|y|]}{\sigma_w}$ or as
$\frac{\sqrt{E[y^2]}}{\sigma_w}$ or as $\frac{E_b}{N_0}$ (a quantity that I like to call the BEND ratio ("bit energy to noise density ratio")), because those definitions make more sense in that particular application, or lead to more memorable formulas. If $\mathsf{SNR}$ is defined to be the BEND ratio, then the bit error probability for coherent demodulation of
binary orthogonal FSK would be $P_e = Q(\sqrt{\mathsf{SNR}}\,)$ (where $Q(\cdot)$ is the complementary standard Gaussian CDF) while for differentially coherent demodulation of differentially encoded PSK, the bit error probability would be $p_e =\frac{1}{2}e^{-\mathsf{SNR}}$, etc.
Skipping over the fine details to the executive summary, large values of  $\mathsf{SNR}$ are better than small values of $\mathsf{SNR}$.

Since the $x[n]$'s are random (variables), so are the $y[n]$'s random
variables since they are sums of the $x[k]$'s. Note that
\begin{align}
y[n] &= x[n]+y[n-1]\\
&= x[n] + \left(x[n-1]+y[n-2]\right)\\
&= x[n]+x[n-1]+\left(x[n-2]+y[n-3]\right)\\
&\ddots\\
&= \sum_{k\leq n} x[k]
\end{align}
Thus, $z[n] = y[n]+w[n]$ is the sum of two independent random variables
and so
\begin{align}\require{cancel}E[(z[n])^2] &= E[(y[n])^2] + E[(w[n])^2]+2E[y[n]]\cancelto{0}{E[w[n]]}\\
&= E[(y[n])^2] + \sigma_w^2\tag{1}\label{1}\end{align}
where $E$ denotes expectation (as the OP suspected) and the
expectation is indeed calculating the mean, not of the random variable $z[n]$ but of the random variable $(z[n])^2$, and $\sigma^2$ is the variance of the AWGN.  So, the definition of $\mathsf{SNR}$ that the
OP is reading about, viz. $\frac{E[y^2]}{\sigma_w^2}$ is the ratio of the two terms in $(1)$. The numerator is the instantaneous power of the
output signal $y$ and the denominator is the instantaneous power of the AWGN.  This definition of $\mathsf{SNR}$ is deemed by the authors of
the paper that the OP is reading to be the most convenient for their
purposes. Others may choose to carry out similar analyses using different definitions of $\mathsf{SNR}$ and their formulas will appear to be different. But, arguments about "My $\mathsf{SNR}$ is larger than your
$\mathsf{SNR}$" are readily resolved by refusing to compare apples with oranges and looking at the important parameters: which system achieves smaller error probability, or needs smaller (peak or average) transmitter power, uses smaller bandwidth, etc. The rest is just semantics.


Are there libraries for extraction of sound wave features?

From the ones I've been using I can recommend:

YAAFE - very pleasant to work with in Python
ESSENTIA - another one I like particularly due to Python integration
aubio
FEAPI
Aquila - friend of mine used it extensively and he likes it a lot

Recently I came across this paper and I believe that this should perfectly answer your question.

  Moffat D. et al - An Evaluation of Audio Feature Extraction
  Toolboxes

They do provide a detailed comparison based on:

Amount of features available
Ease of usage
Output format of the features
Computational efficiency

They tested following libraries:

Aubio
Essentia
jAudio
Librosa
LibXtract
Marsyas
Meyda
MIR Toolbox
Timbre Toolbox
YAAFE


Absolute maximum number of possible PN sequences for given length

No, that isn't possible. Think about it in the geometric sense, where each bit in the sequence is a dimension. In that context, you can recast your question as:

  In an $n$-dimensional space, are there any sets of $M$ vectors (where $M > n$) that are mutually orthogonal?

For example, in three-dimensional space, can you construct a set of four vectors that are mutually orthogonal to one another? You can't. Stated in linear algebraic terms, you can only have as many basis vectors as the dimensionality of the space.


Relation between constellation SNR and baseband SNR

I think that the question by OP xvan has several misconceptions in it and these propagate into the answer that he has provided for his own question.

Define
\begin{align}\phi_1(t) &= \sqrt{\frac 2T}\cos(2\pi f_ct),~ 0 \leq t \leq T,\\
\phi_2(t) &= -\sqrt{\frac 2T}\sin(2\pi f_ct),~ 0 \leq t \leq T,
\end{align}
where $f_c$ is the carrier frequency are two orthogonal unit energy RF signals and thus $$s(t) = s_{I}\phi_1(t)+s_{Q}\phi_1(t)\tag{1}$$ is a RF QPSK signal with rectangular baseband pulses (unless, of course, $s_I$ and $s_Q$ are themselves some suitably chosen pulse shapes (e.g. raised-cosine) whose dependence on $t$ has been suppressed by the OP.  (Note that I have inverted $\phi_2(t)$ for reasons that will become apparent soon).  $s(t)$ as defined in $(1)$ is a real-valued signal and one does not add complex Gaussian noise to it; one must add (real-valued continuous-time) white Gaussian noise to $s(t)$ to get $x(t)$.  It is also helpful to assume that $f_cT$ is an integer, and if not, at least that $f_T \gg 1$.
The complex-valued baseband symbol $[s_1, s_2]$ (or $\mathbf s = s_1+js_2$ if one likes) takes on values $\left[\pm \sqrt{\mathcal E_b}, \pm \sqrt{\mathcal E_b}\right]$ where $\mathcal E_b$ is the signal energy per bit. These are the constellation points in the complex baseband equivalent of the RF signal $(1)$. Thus,
\begin{align}
s(t) &= \pm\sqrt{\frac{2\mathcal E_b}{T}}\cos(2\pi f_ct)
\mp \sqrt{\frac{2\mathcal E_b}{T}}\cos(2\pi f_ct)\\
&= \Re\left(\mathbf s \sqrt{\frac 2T}\exp(j2\pi f_c t)\right) \tag{2}
\end{align}
showing that the QPSK signal is the real part of the result of (amplitude-)modulating the complex baseband signal $\mathbf s$ onto the complex carrier $\sqrt{\frac 2T}\exp(j2\pi f_c t)$.
Note also that if I had stuck with the OP's definition of $\phi_2(t)$, $(2)$ would have been
$$s(t) = \Re\left(\mathbf s \sqrt{\frac 2T}\exp(-j2\pi f_c t)\right)$$
which would have been a disaster since I don't believe in negative frequencies on Mondays, Wednesdays and Fridays.
To the signal in $(1)$, we add continuous-time white Gaussian noise $W(t)$ which is a hypothetical zero-mean process with autocorrelation function $R_W(t) = \frac{N_0}{2}\delta(t)$ where $\delta(t)$ is the Dirac delta, that is, the power spectral density (PSD) $S_W(f)$ has value $\frac{N_0}{2}$ for all $f, -\infty < f < \infty$. Hypothetical because the noise has infinite power. More practically, the result of passing $W(t)$ through a BIBO filter with transfer function $H(f)$ results in a zero-mean WSS Gaussian process with PSD $\frac{N_0}{2}|H(f)|^2$ while
$\int_{0}^{T} \phi_i(t)W(t) \,\mathrm dt$ is a zero-mean Gaussian random variable with variance
\begin{align}
\sigma^2 &= E\left[\int_{0}^{T} \phi_i(t)W(t) \,\mathrm dt\int_{0}^{T} \phi_i(u)W(u) \,\mathrm du\right]\\
&= E\left[\int_{0}^{T} \int_{0}^{T}\phi_i(t)\phi_i(u)W(t) W(u) \,\mathrm dt \,\mathrm du\right]\\
&= \int_{0}^{T} \int_{0}^{T}\phi_i(t)\phi_i(u)E[W(t)W(u)] \,\mathrm dt \,\mathrm du\\
&= \int_{0}^{T} \int_{0}^{T}\phi_i(t)\phi_i(u)\frac{N_0}{2}\delta(t-u) \,\mathrm dt \,\mathrm du\\
&= \int_{0}^{T} \phi_i(u)\left[\int_{0}^{T}\phi_i(t)\frac{N_0}{2}\delta(t-u) \,\mathrm dt\right] \,\mathrm du\\
&= \frac{N_0}{2}\int_{0}^{T} \phi_i^2(u)\,\mathrm du\\
&= \frac{N_0}{2}
\end{align}
since $\phi_1(t)$ and $\phi_2(t)$ are unit-energy signals. A calculation similar to the one above shows that $\int_{0}^{T} \phi_1(t)W(t) \,\mathrm dt$ and $\int_{0}^{T} \phi_2(t)W(t) \,\mathrm dt$ are uncorrelated random variables and hence independent random variables since they are jointly Gaussian random variables.  In short, for $i=1,2$,
\begin{align}
\hat{s}_i &= \int_{0}^{T} x(t)\phi_i(t) \,\mathrm dt\\
&= \int_{0}^{T} [s(t) + W(t)] \phi_i(t) \,\mathrm dt\\
&= s_i + W_i
\end{align}
where $s_i$ has value $\pm \sqrt{\mathcal E_b}$ while $W_1$ and $W_2$ are independent Gaussian random variables with identical variance $\frac{N_0}{2}$.

Now, the OP wants to sample $x(t)$ a large number of times (perhaps multiple samples per period of the carrier) and use the samples presumably to simulate the operations of the correlators in the receivers to make sure that the correlators are doing what the math says that they are doing. This is a whole different kettle of fish since one cannot sample white noise in any meaningful sense and what the variance of the samples turns out to be depends on the details of the actual sampling device and not any model of the device that one might devise in MATLAB. If the "sampler" effectively grabs a very very short segment of $x(t)$ and returns the average value of $x(t)$ over that segment, then model that, and not that the sampler puts out the instantaneous value of $x(t)$.

But, let's take a more realistic/charitable view and assume that the noise process is band-limited white noise whose PSD has value $\frac{N_0}{2}$ only for $|f| \in \left[f_c-\frac 1T, f_c+\frac 1T\right]$ which is the nominal null-to-null passband for the QPSK signal. The noise autocorrelation function $R_W(t)$ is thus $\frac{N_0}{2}\operatorname{sinc}\left(\frac tT\right)\cos(2\pi f_ct)$ from which the most important take-away is that samples from $W(t)$ that are spaced $T$ seconds apart are independent, but samples spaced less than $\frac{1}{2f_c}$ seconds apart are very definitely correlated.  The samples are indeed Gaussian random variables and thus the Riemann sums that be used are indeed Gaussian random variables, but we can't just add up the variances of the variables in the sum; they are correlated random variables and so we must take the covariances into account as well. This messes up the calculations and getting it all straightened out is beyond what one should expect in an answer in a forum such as dsp.SE. What the OP has written in his own answer is not quite right, but is not correctable in a few words.


How to obtain the exact value of wavelength from a 2D FFT amplitude vs wavenumber plot like it is obtainable from 1D FFT amplitude vs wavenumber plot?

First, you want to find the peaks in this 2D plot. In your case, you know that there will be 2 peaks in your FFT (because you added two sines). Your plot shows 4 peaks because the FFT magnitude is symmetric, since your signal is real. You can discard negative wave numbers and look for peaks among the positive wave numbers.
The wave numbers corresponding to those peaks belong to your signal. There will be two wave numbers corresponding to each peak, $k_x$ and $k_y$. We know that $k = \frac{2\pi f}{c}$, where $k = \sqrt{k_x^2 + k_y^2}$ is the wave-number, $f$ is the frequency in Hz and $c$ is the speed of sound in the medium. Now, using $\lambda = \frac{c}{f}$, we can deduce that $\lambda = \frac{2 \pi}{k}$.


For the case of linear modulation on the AWGN channel with equiprobable symbols (a very common case) the optimum approach is to truly use a filter that is matched to the symbol waveform, i.e.:

$$
q(x) = p(x)
$$

Using a matched filter provides the optimum signal-to-noise ratio at the filter output at each decision instant. This is easy to see when you remember that a matched filter acts like a sliding cross-correlator between its input signal and the expected symbol waveform, correlating the two at all possible lags. At optimum decision instants, the filter's impulse response (typically scaled to have unit energy) lines up exactly with a transmitted symbol, analogous to a zero-lag condition on the cross-correlation operation. At this time value, the filter's output is equal to the amount of energy in the received symbol, scaled by a data-dependent factor (e.g. for BPSK, the matched filter would output $ E_s $ or $ -E_s $), plus a noise term.

The noise energy at the filter output during the sampling instant is not dependent upon the time-domain shape of the filter's impulse response, only the impulse response's total energy (as noted previously, typically unity). Therefore, the signal to noise ratio is maximized by maximizing the amount of signal energy in the filter output at the sampling instant. By choosing the receiver filter to be matched to the symbol shape, we have done so, as the symbol waveform has maximum correlation with a filter impulse response that has an identical shape. Thus, the matched filter provides maximum SNR, for the AWGN channel case.

With that bout of hand-waving out of the way (you can definitely get at it with more mathematical rigor, but I'm an engineer and this is a free service; if you want to dig into the details, check any digital communication theory text), you might be thinking that I forgot that you asked about the non-ideal, ISI case. Fear not, for I assert that if you know the transmitted pulse shape, the matched filter is still the optimum choice for the AWGN channel.

The key: if you know the responses of the pulse-shaping and receiver detection filters $p(x)$ and $q(x)$, and the last "few" transmitted symbols, you can calculate what the ISI induced by those previous symbols would be and account for it accordingly; it is a deterministic quantity. The amount of symbol history that you require is related to the amount of ISI you have, i.e. how many symbol periods the cascaded filter response smears across.

Of course, you typically don't know with certainty what the previous few symbols were; if you did, then you might be at a high enough SNR that your ISI can be neglected. In the more interesting case, you can't make that assumption. Instead, a maximum-likelihood sequence detection approach is employed using the Viterbi algorithm. This process is referred to as Viterbi equalization, because in this model you treat the ISI induced by the pulse shape like a soft-valued convolutional code that is applied to your transmit waveform. The time duration of the ISI in the Viterbi equalizer defines the required number of algorithm states, similar to the constraint length in a convolutional code.

This approach is often used in systems that have the non-optimum pulse shape that you noted; one notable example is GSM (which uses a Gaussian pulse shape that extends across multiple symbol intervals). One great reference on this topic was published by Sklar in 2003:

B. Sklar, “How I learned to love the trellis“, IEEE Signal Processing Magazine, pp. 87-102, May, 2003


Recover curves from noisy collection of points

Your problem is very similar to the cell-tracking problem. This has been solved quite well using tracklet based approaches, where each object (bubble, cell etc) is modeled as a Maximum-a-Posteriori (MAP) estimation. The solution to which is through linear programming.
You can get a good idea of the approaches and a web-based implementation here


Here is a rough cut.  You could definitely get more elaborate on your peak testing, e.g. comparing it to neighbor peaks.  This simply makes a metric for how far the peak "sticks out" and uses a cutoff value on that and the height to select a candidate peak.

In the case of two or more  isolated peaks in your field of zeroes, it would select all of them.  Not sure that is what you want, but this works for your given example.

Once your peaks are selected, then you can zero all other data values.

Ced

        thePeak_Spot       = np.zeros(100, dtype=int )
        thePeak_Height     = np.zeros(100)
        thePeak_Protrusion = np.zeros(100)

        thePeakCount = 0

        for n in range( 1, len( data ) - 1 ) :

            theLeftDrop  = data[n] - data[n-1]
            theRightDrop = data[n] - data[n+1]

            if theLeftDrop > 0 and theRightDrop > 0 :
                thePeak_Spot[thePeakCount] = n
                thePeak_Height[thePeakCount] = data[n]
                thePeak_Protrusion[thePeakCount] = (theLeftDrop + theRightDrop) / ( 2.0 * data[n] )
                thePeakCount += 1

        for p in range( thePeakCount ) :
            if thePeak_Protrusion[p] > 0.5 and thePeak_Height[p] > 0.1 :
               print p, thePeak_Spot[p], thePeak_Height[p], thePeak_Protrusion[p]

        plt.plot(data)
        plt.show()

output:

4 13 0.497254209953 0.695280339927
5 58 0.182008917152 1.0
6 122 0.288306083312 0.900874833979


Why do AAC files seem to have an inexact bitrate?

Both AAC and MP3 are not always encoded with constant bit rate but they support variable bit rate encoding too.


Differences between matlab's 'impulse' and 'lsim' command; create impulse response

Indeed the lsim function is returning scaled values. You might also notice that it is connected to the gain $k$ of your filter:

[z,p,k] = tf2zpk(sysD.num{1}, sysD.den{1});

That will return:

k=0.002

You can scale it back by the sampling period as you noticed.


I think you have a difficult problem here. Smoothing you signal (for example, with a moving average) probably won't help you. I guess you could try passing your signal through a narrowband notch filter whose notch center frequency is the frequency of your unwanted pulses. MBaz is correct. To think about an appropriate filter means you'd have to know what your signal looks like if it were not contaminated by those unwanted pulses. For example, does your uncontaminated signal look like the following?

I realize I'm not answering your question in any meaningful way, but I was thinking: I assume your noisy signal is the output of an A/D converter sampling an analog voltage signal that's badly contaminated by some sort of electrical noise, some electrical interference. Instead of eliminating the noise in the digital samples I think you should do everything you possibly can to eliminate the electrical noise in your original analog signal.

Dec. 12, 2015 Update based on Vlad's Comment:
@Vladimir: Ah, interesting! What you wrote makes much sense. Those pulses are not noise contamination, they are part of your signal. To attenuate those pulses, three things off the top of my head:

[1] I'd like to see the result of applying your signal to a moving average process.

[2] Although I have no meaningful experience with such a filter, I wonder what would be the result of applying your signal to a median filter.

[3] There's a scheme in the "DSP Tricks" chapter of my DSP book describing how to attenuate very narrow impulsive noise in a signal. I don't know if that scheme would be useful with regard to your pulses but who knows? If you'd like me to send you information on this scheme, send me a private e-mail at: R_dot_Lyons_@_ieee_dot_org.


How to find period/envelope that is larger than sample size?

If you can guess the number of Fourier coefficients, then a decomposition method such as MUSIC or ESPIRIT might provide a parametric estimate or eigenvector decomposition, given that guess.


Need to calculate $\frac{E\left[|x(n)|^4\right]}{E\left[x^2(n)\right]}$

The most efficient way that I can think of to do it would be as follows:

For each input sample $x[n]$, calculate $|x[n]|^2$. This is more efficient than calculating $|x[n]|$, because you don't have to do the square root:

$$
|x[n]| = \sqrt{\text{Re}\{x[n]\}^2+\text{Im}\{x[n]\}^2} \\
|x[n]|^2 = \text{Re}\{x[n]\}^2+\text{Im}\{x[n]\}^2
$$

Average the instantaneous power estimate $|x[n]|^2$ as needed for your application (you mentioned exponential averaging).
Note that you can also get $|x[n]|^4$ efficiently using this method; you just square the measurement of $|x[n]|^2$ that you made above, then average again as desired.


I think your recalculation of phis in the function is not correct.

If I do this in scilab:

//25631
N = 100;
phis1=162*%pi/180;//initial phase of signal

[data1,phis2] = rc(N,phis1);
[data2,phis3] = rc(N,phis2);
[data3,phis4] = rc(N,phis3);

clf
plot([data1 data2 data3])

where rc is

////////////////////////////////////////////////////////////////////////
//  generates a sinusoid of dwLength in dwData
//      and maintains phis value updated
////////////////////////////////////////////////////////////////////////
function [dwDataOut,phis] = rc(dwLength, phis)
omegaf=50/dwLength/%pi*90/100;

dwDataOut = zeros(1,dwLength);

for dwi=1:dwLength,
    dwDataOut(dwi) = 4450*sin(dwi*omegaf+phis)+200;
end

//update wave phase
//WAS: phis = phis - %pi*(1-50/dwLength)+%pi/2;
phis = dwLength*omegaf+phis;
while(phis<0)
    phis = phis + 2*%pi;
end
endfunction

Note specifically the line phis = dwLength*omegaf+phis;

Then the plot is as below, which seems to be what you're after.


How to find the impulse response of this system with complex roots?

It turns out that those terms with complex conjugate scalers and exponentials can be combined.  Rather than do it to your example, I would rather do it in the Z-domain.  Let

$$\begin{align}
H(z) &= \frac{A_1}{1-p_1\,z^{-1}} + \frac{A_2}{1-p_2\,z^{-1}} \\
\\
&= \frac{A_1(1-p_2\,z^{-1})+A_2(1-p_1\,z^{-1})}{(1-p_1\,z^{-1})(1-p_2\,z^{-1})}  \\
\\
&= \frac{(A_1-A_1 p_2\,z^{-1})+(A_2-A_2 p_1\,z^{-1})}{1-(p_1+p_2)\,z^{-1} + p_1 p_2\, z^{-2}}  \\
\\
&= \frac{(A_1+A_2) - (A_1 p_2 + A_2 p_1)z^{-1}}{1-(p_1+p_2)\,z^{-1} + p_1 p_2\, z^{-2}} \\
\\
\end{align}$$

Now it turns out that there are two ways for the impulse response, $h[n] = \mathcal{Z}^{-1} \{H(z)\}$, to be total real (with no imaginary component).  One way is for both poles $p_1,p_2$ to be real and that will also require the two coefficients $A_1,A_2$ to also be real.

$$\begin{align}
h[n] &= \mathcal{Z}^{-1}\{H(z)\} \\
\\
 &= \mathcal{Z}^{-1}\left\{ \frac{A_1}{1-p_1\,z^{-1}} \right\} + \mathcal{Z}^{-1}\left\{ \frac{A_2}{1-p_2\,z^{-1}} \right\} \\
\\
&= A_1 p_1^n u[n] + A_2 p_2^n u[n] \\
\end{align}$$

The other way for $h[n]$ to be real is for both poles $p_1,p_2$ to be complex-conjugate of each other and that will also require the two coefficients $A_1,A_2$ to also be complex-conjugate.

$$\begin{align}
h[n] &= \mathcal{Z}^{-1}\{H(z)\} \\
\\
 &= \mathcal{Z}^{-1}\left\{ \frac{(A_1+A_2) - (A_1 p_2 + A_2 p_1)z^{-1}}{1-(p_1+p_2)\,z^{-1} + p_1 p_2\, z^{-2}} \right\}  \\
\end{align}$$

(i need to return to this later.)


This sounds counter intuitive to many, but as long as the difference operator and the smoothing kernel are linear and space-invariant, they can be applied in any order, and thus are often combined in a single convolution operator (for more computational efficiency), for the same.
For some intuition, consider that, either for the linear smoothing and the derivative, a noisy pixel  is replaced by a linear weighed combination of other pixels. And a linearly-weighed combination of another linearly-weighed combination remains a linear combination, the same due to the commutativity and the associativity properties of the addition and the multiplication.
This is illustrated below for the first derivative only, but you can follow the same argument by differentiating again.

What you describe is commonly known as the LoG, Laplacian of Gaussian.
Note that when one operation becomes nonlinear, this is not true anymore.


MLE parameter estimation -- confusion regarding some terms in the pdf of complex normal r.v (Part 2)

Hi: It's worded more clearly now in that you're estimating A and there's only one RV which makes more sense. Consider the first link I sent in the other message and go to where it says the "likelihood is". Your likelihood is quite similar except, since your variance is multiplied by 2, this causes 2 things to happen in terms of how would you change the likelihood in that document accordingly for your problem.

A) The $2 \pi \sigma^2$ that is getting raised to the $-\frac{n}{2}$ instead becomes  $4 \pi \sigma^2_{v}$

B) The denominator of the exp instead of  being $2 \sigma^2$  becomes $4 \sigma^2_{v}$.

Make those 2 changes to that likelihood and you will have your expression. I get confused looking at your 2) because my memory of 1/2's and the squares ( in your case H's ) is blurry and it can be confusing. So better to rely on document where likelihood is correct so you can't go wrong. Keep in mind that this is under the assumption that complex RV's don't change much compared to real ones.


I do not know what convolve.convolve does but the output of signal.convolve is the linear convolution (as opposed to circular convolution) of the two sequences. And since you supplied 'same', the output has the same size as the first input vector, and the samples are taken from the center of the complete output.

The circular convolution is indeed

34   32   34   40

and the linear convolution is

5   14   26   40   29   18    8

from which you obtained samples 2 to 5 (by using the mode 'same').


how to detect center of a blurry circle with opencv

You'll need to improve the contrast of the image first, then filter it strongly to remove the noise. since the circles are 'thick' (blurry), you can filter quite a bit without destroying the circles structure.

I would then apply some edge detection algorithm to get a binary edge image that can be processed by the Circular Hough Transform.

I get the following edge image from your image:

Using the following MATLAB commands:

 % x is the input grayscale image. First we adaptively improve the contrast over the image
 y= adapthisteq(x);

 % next we use the Canny edge detector with a strong Gaussian lowpass filter
 ee=edge(y, 'canny', [], 5);


Clarification Regarding Changing the Variable

The implication seems a little bit wrong, but that is okay.
When dealing with the functions in the frequency domain, let’s discuss it over the function $X(j\omega)$, when shifting or scaling the function on the frequency axis, the operations are done on the angular frequency $\omega$ or the frequency $f$.
On the frequency axis, there is no imaginary number, that is, the frequency axis is full of real numbers. Therefore, calling a point on the frequency axis as e.g. $j2$ Hz is incorrect.
When I do calculations on the paper, I prefer using $\omega$ instead of $j\omega$ when defining the variable of the Fourier transform of a function e.g. $X(\omega)$. However, in the case of transferring a Fourier transform into the Laplace transform under specific conditions, $j\omega$ notation aids in preventing confusion of where to put the $s$ variable as in this process, we do assignment $s = j\omega$.
I've changed the given function a little bit i.e. I've taken the absolute value of the angular frequency $\omega$ in the power of the exponential in order to illustrate what happens better.
Okay, along with these knowledge, let’s focus on your function. The amplitude spectrum of the function $X(j\omega)$ is given below:

As it is noted, the mathematical explanation of the function is;
$$X(j\omega)=\begin{cases}e^{-2|\omega|},&\text{if $|\omega|<2\pi$}\\0,&\text{otherwise}\end{cases}$$
So, let’s shift the function $X(j\omega)$ by $\Omega$ to the right on the frequency axis as I’ve assumed $\Omega>0$. This shifting process, which is a frequency shifting, is done in the same way as time shifting. The resultant shifted function $X(j(\omega - \Omega))$ will be in the figure below:

So, the resultant mathematical expression becomes;
$$X(j(\omega - \Omega))=\begin{cases}e^{-2(|{\omega - \Omega}|)},&\text{if $|{\omega - \Omega}|<2\pi$}\\0,&\text{otherwise}\end{cases}$$
The expression $|{\omega - \Omega}|<2\pi$ can be represented comprehensively as ${{-2\pi}+{\Omega}}<{\omega}<{{2\pi}+{\Omega}}$.


Need for Oversampling in OFDM and symbol time determination

As the questioner refers to my code at github in the comments, I'd like to clarify.

There is no need to apply oversampling in OFDM, nor is oversampling required to recover input signals. However, in practical implementations you have to consider some additional constraints regarding the signal's bandwidth. The WARP v2 radio board is clocked with 40 MHz, i.e., the Nyquist frequency is 20 MHz. However, the WARP v2 architecture limits the generated baseband signal's bandwidth to ~10 MHz by TX and RX lowpass filters. This is a meaningful measure to avoid unpleasant aliasing effects. Hence, the OFDM signal's bandwidth should be limited to 10 MHz in the baseband. You can easily achieve this by padding all subcarriers above 10 MHz with zeros. This will work perfectly well.
An alternative (but probably rather uncommon) way to reduce the bandwidth is to apply 2x oversampling while using the majority of all subcarriers. I implemented oversampling in my generic (not-OFDM-related) transmission procedure and finally enabled it in my OFDM code as I could achieve a slightly higher SNR in our office environment with infrequent interference and some narrowband interference. However, you can safely disable oversampling in my code and apply zero-padding instead.
An OFDM symbol is the time-domain representation of the symbols concurrently transmitted over all subcarriers. It consists of just as many samples as there are OFDM subcarriers since it is calculated as the IDFT of all subcarriers' symbols. When adding the cyclic prefix, you add $N_G$ samples. The duration in seconds is the number of samples times the time per sample. When oversampling, the symbol duration extends according to the oversampling factor.

Kind regards,
Robin


Can all systems be described by differential/difference equations?

No the converse is not true in general.

Take for example the discrete-time ideal lowpass filter with impulse response

$$ h[n] = \frac{ \sin( \omega_c n) }{ \pi n } ~~~,~~~-\infty < n <\infty$$

which describes an LTI system but it does not correspond to a difference equation of any kind. Indeed $h[n]$ is derived based on the inverse discrete-time Fourier transform of the filter frequency response;

$$ H(\omega) = \begin{cases} {1 ~~~,~~~ |\omega| < \omega_c \\ 0 ~~~,~~~ \text{ otherwise } }\end{cases} $$

$$h[n] = \frac{1}{2\pi} \int_{-\pi}^{\pi} H(\omega) e^{j\omega n} d\omega $$


Is there any good book regarding 3D reconstruction from 2D images? I prefer to work in MATLAB

There are a lot of good books. Multiple View Geometry in Computer Vision,

"Computer Vision: Algorithms and Applications", but anyway, replication of the theory is the necessary stage to me.

You can start from the finishing course "3D Computer Vision" from the CTU. They published lectures and labs - in MatLAB.

P.S. Here is a ready implementation http://www.fit3d.info/.


The circular convolution is actually mainly an artefact from the Fast Convolution algorithm, that is calculating the convolution via the FFT. I am not aware of any use of the circular convolution at all. Now, to efficiently calculate the convolution between x and h, you will often use the FFT, as a multiplication in the frequency domain equals a convolution in the time domain (and inversely). The algorithm is then:

conv(x,h) = ifft( fft(x) .* fft(h) )

where .* denotes an element-wise multiplication. You will have to use a zero-padding so both vectors have the same length and you can do the multiplication. This way to calculate the convolution is faster than the direct way in the time domain. However, as fft(x), fft(h) and thus also ifft(...) have the same length (here, the length of x) you cannot express the "real", linear convolution, which would be a longer vector. What happens is, that the result of the linear convolution gets "wrapped around" as you indicate in your drawings.

Let me show you that using your example. The linear convolution is (without proof)

conv(x,h) = [2, 5, 8, 8, 5, 4, 4]

Now as in the circular convolution, the vectors are only 5 long, the last 2 entries of the result will be added in the front, so

circ_conv(x,h) = [2+4, 5+4, 8, 8, 5] = [6, 9, 8, 8, 5]

is the circular convolution. The easiest way (imho) is to first calculate the linear convolution and then wrap around that result to achieve the circular convolution.

If you want to calculate the linear convolution with the FFT, you only have to zero-pad both vectors such that the linear convolution fits into these vectors, i.e. to a length of 7 in your example.