transient and steady-state response of first order system Considering this general 1st order transfer function

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

How to find (analytically) the transient and steady-state responses?

With steady-state response I mean response to a sinusoid. I'm particularly interested in the case when $a=\pm1$.

Conceptually I have some difficulties grasping how to deal with the fact that the z-transform for a causal sinusoid only exists for $|z|>1$ (explained here) and that for the case when $a=\pm1$ the z-transform does not converge for $|z|=1$. So the 'straight-forward' procedure of multiplying $H(z)$ with the z-transform of a causal sine and then inverse z-transform the result seems to not be possible/valid. I'm confused.

Edit:

simple example

fs = 48000;
f0 = 1000;
T = 3/f0;
N = round(T*fs);
t = linspace(0,T,N);
x = sin(2*pi*f0*t);

bd = [1 0];
ad = [1 -1];
yi = filter(bd, ad, x);
figure;
plot(t,x,'b'); hold on;
plot(t,yi,'r')

The region of convergence (ROC) of the $\mathcal{Z}$-transforms of a step-modulated sinusoid

$$x[n]=\sin(\omega_0n)u[n]\tag{1}$$

is the region $|z|>1$. The ROC of the transfer function

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

also equals $|z|>1$ (assuming a causal system). Consequently, multiplying the two transforms doesn't pose any problem, and the result also converges for $|z|>1$.

As an example, let's compute the response $y[n]$ of the system $(2)$ to the input $(1)$. For the sake of simplicity we assume $b_0=1$ and $b_1=0$ in $(2)$. First, we compute the response to

$$\tilde{x}[n]=e^{j\omega_0n}u[n]\tag{3}$$

The $\mathcal{Z}$-transform of $(3)$ is

$$\tilde{X}(z)=\frac{1}{1-e^{j\omega_0}z^{-1}}\tag{4}$$

The $\mathcal{Z}$-transform of the response $\tilde{y}[n]$ is given by

$$\begin{align}\tilde{Y}(z)&=\tilde{X}(z)H(z)\\&=\frac{1}{1-e^{j\omega_0}z^{-1}}\cdot \frac{1}{1-z^{-1}}\\&=\frac{A}{1-e^{j\omega_0}z^{-1}}+\frac{A^*}{1-z^{-1}}\tag{5}\end{align}$$

with

$$A=\frac{1}{1-e^{-j\omega_0}}=H(e^{j\omega_0})\tag{6}$$

The inverse $\mathcal{Z}$-transform of $(5)$ is

$$\tilde{y}[n]=H(e^{j\omega_0})e^{j\omega_0n}u[n]+H^*(e^{j\omega_0})u[n]\tag{7}$$

The response to the step-modulated sinusoid $(1)$ is easily obtained from $(7)$ by taking its imaginary part:

$$\begin{align}y[n]&=\textrm{Im}\big\{\tilde{y}[n]\big\}\\&=\big|H(e^{j\omega_0})\big|\sin\big(\omega_0n+\arg\left\{H(e^{j\omega_0})\right\}\big)u[n]\\&\qquad -\textrm{Im}\big\{H(e^{j\omega_0})\big\}u[n]\tag{8}\end{align}$$

The first term in $(8)$ is the steady-state response, and the second term is the transient response, which doesn't decay because of the system's pole at $z=1$.

So the output has a DC component due the imaginary part of $H(e^{j\omega_0})$. Note that this is no contradiction with the linearity of the system, because the input signal is not a single spectral line but it has a continuous spectrum extending down to DC due to the sinusoid being switched on at $n=0$. This DC value of the input's Fourier transform triggers the system's eigenfrequency, which is a DC component. Note that this always happens with LTI systems: the output may contain oscillations at frequencies that are only determined by the system, not by the input signal. However, unlike in the given example, usually these transients decay because most of the time we consider asymptotically stable systems.

I slightly modified your Matlab/Octave script and added the analytical result $(8)$ for comparison. The analytical result is identical (up to numerical accuracy) to the result obtained by filtering:

fs = 48000;
f0 = 1000;
w0 = 2*pi*f0/fs;
N = 200;
n = 0:N-1;
x = sin(w0*n);

bd = [1 0];
ad = [1 -1];
yi = filter(bd, ad, x);
figure;
plot(n,x,'b'); hold on;
plot(n,yi,'r')

% analytical computation
A = 1 / ( 1 - exp( -1i*w0 ) );
y2 = abs(A) * sin( w0*n + angle(A) ) - imag(A);

plot(n,y2,'k.'), hold off
legend('input','response by filtering','response by analytical computation')


What is the difference between multiplying a delta and a step versus convolving a delta and a step? Seems both will produce another step. there is no difference? Thanks

First of all you need to see whether you are performing these operations for a continous time signal or discrete time signal.

Sampling theorem says that multiplication of a signal $x(t).\delta(t)$=$x(0).\delta(t)$ provided $x(t)$ is continous at $t=0$. But here in your question $x(t)$ is a unit step function which is not continous at $t=0$. Hence the multiplication of $\delta(t).u(t)$ is not defined.

In the discrete case the same property is $\delta[n].x[n]$=$x[0].\delta[n]$ and there is no question of continuity as signal is discrete. In this case your question can be written as $u[n].\delta[n]$=$u[0].\delta[n]$=$\delta[n]$ as $u[0]$ is equal to 1.

Now the convolution property says that $x(t)*\delta(t-t_0)$=$x(t-t_0)$ where $*$ is denoted here as the convolution operator. In your question this would be $u(t)*\delta(t)$=$u(t)$ which is nothing but unit step function for continous time.

The same is true for discrete signals i.e.$x[n]*\delta[n-n_0]$=$x[n-n_0]$. So for discrete case your question would be $u[n]*\delta[n]=u[n]$ which is a unit step function for the discrete case.

You can look at comments of Matt below my answer for better clarity.


Constructing a tone through additive synthesis This question follows from: Synthesizing harmonic tones with Karplus Strong

I am attempting to generate a full soundbank of 88 notes.

As far as I can see my three options are:

Using samples
Using something like Karplus Strong (as linked above)
Synthesising each harmonic and combining them

(1) does not give much flexibility and (2) is difficult: to gain any kind of control requires a sophisticated level of understanding of control systems, and I don't have any training. I've been looking through an online MIT lecture series, and I'm daunted. This would take me months or even years to understand.

One disadvantage of (3) is that it may not handle real-time composition as well as (2); it may be heavier to process. For me that's not a problem, I don't need real-time.

I think (3) should be much easier than (2). It should be possible to construct a formula for the amplitude A(t) of a particular harmonic.  But I imagine there must be some science of shaping these harmonics; I'm sure many people have experimented with this, and there must be publicly available literature.

And I'm guessing there are some techniques for improving the quality of the sound, for example maybe giving some A * exp( kt ) sin( wt ) wobble to each harmonic, and maybe a note that is struck louder could have a correspondingly higher value for A.

etc.

So my question is: what is the parameter space? i.e.  what are the best things to tweak?

Links to relevant resources are gratefully received!

π

PS please note that I'm not after any particular timbre, certainly I'm not try to emulate a piano, I'm simply after a harmonic tone that is pleasing to the ear.

You say:

  I'm not after any particular timbre, certainly I'm not try to emulate a
       piano, I'm simply after a harmonic tone that is pleasing to the ear.

but I think something that is pleasing to the ear is probably harder to quantify than emulate a piano!

Bill Sethares has done lots of good work (both theoretically in signal processing and practically as a musician).  Bill wrote a book specifically about these ideas.

You can hear an example of his work, called C to shining C here. There is a discussion about it on Wikipedia.

The central idea is that consonance and dissonance are related to psychoacoustics: how the brain perceives sound.  As Bill says:

  For instance, the perception of "timbre" is closely related to (but also quite distinct from) the physical notion of the spectrum of a sound. Similarly, the perception of "in-tuneness" parallels the measurable idea of sensory consonance. The key idea is that consonance and dissonance are not inherent qualities of intervals, but are dependent on the spectrum, timbre, or tonal quality of the sound. To demonstrate this, the first audio track on the accompanying CD plays an example where the octave has been made dissonant by devious choice of timbre, even though other, non-octave intervals remain consonant. In fact, almost any interval can be made dissonant or consonant by proper sculpting of the timbre.

So, now that I've said it's too hard to answer your question, I'll try to give some pointers.

The main parameters to look at are the amplitudes ($A_k$), decay parameter ($\alpha_k$), (relative) phases ($\phi_k$) and frequency ($\omega_k$) relationships between your tones:

$$
s(t) = \sum_{k=0}^{K-1} A_k e^{\alpha_k t} \sin( \omega_k t + \phi_k)
$$

For a harmonic sound, you probably want to select $\omega_k = (k+1) \omega$ where $\omega$ is your fundamental frequency.

For "brighter" sounds, you want $A_k$ for $k>1$ to be larger (i.e. the higher harmonics have greater amplitude).

You can do flanging or phasing by modulating the phase:
$$
\phi_k(t) = \Phi_k \sin(\omega_\phi t)
$$
where $\omega_\phi$ is much lower frequency than $\omega_0$.


How to separate a low amplitude signal of same frequency but different phase I have a signal consisting of two sine waves of same frequency but have different phase (and also a lot of HF noise). Their amplitudes differ by about $120\textrm{ dB}$. How to I go about separating them? SNR is about $130\textrm{ db}$ but its frequency is far higher than my signal, i.e. $f(\textrm{noise}) \gg  f(\textrm{signal})$.

Any special technique? Also I am really just a beginner, so extra points for those who give a detailed explanation.

Two exact sine waves of the same frequency but different phases sum to another exact sine wave.  That resultant sine wave can be decomposed into an infinite numbers of pairs (or any other greater number) of sine waves of the same frequency, not just the two original ones.  Thus more information is required to reduce the possible solution space below infinite.


Is there an equation for efficiency of a 2-radix FFT as you raise $k$? I know the efficiency for an $N$ point 2-radix FFT is $N\log_2(N)$ but assuming $k\leq N$, what if you were looking for the efficiency of calculating $k$ positions of the FFT? Would the efficiency be $k\log_2(N)$?

No.

The FFT is pretty much an "all or nothing" algorithm. There is almost no efficiency gain in reducing the number of desired outputs.

For $k \ll N$ the Goertzel algorithm is typically the better choice. See https://en.wikipedia.org/wiki/Goertzel_algorithm. The break even point is quite small. Maybe $k = 4\cdot log_2(N)$ or thereabouts.


Feedforward disturbance signal compensation The reference Feedforward control concepts through Interactive Tools says, to minimise the effect of the disturbance $d$ on the process output $y$, a compensator $G_{\text{ff}}=P_3/P_1$ shall be used. From a theory point of view, no problem. Even a few simple examples behave as expected.

Encouraged by that, I decided to implement three-phase inverter control in $dq$ reference frame. Inverter model is represented by LC filter equations:
$$
\text{L} \frac{\operatorname{d}}{\operatorname{d}t}
\begin{bmatrix}
i_\text{d} \\
i_\text{q} \\
\end{bmatrix} =
\begin{bmatrix}
v_\text{dn} \\
v_\text{qn} \\
\end{bmatrix} -
\begin{bmatrix}
v_{\text{d load}} \\
v_{\text{q load}} \\
\end{bmatrix} + \text{L} \omega
\begin{bmatrix}
i_\text{q} \\
-i_\text{d} \\
\end{bmatrix}
$$
$$
\text{C} \frac{\operatorname{d}}{\operatorname{d}t}
\begin{bmatrix}
v_\text{d load} \\
v_\text{q load} \\
\end{bmatrix} =
\begin{bmatrix}
i_\text{d} \\
i_\text{q} \\
\end{bmatrix} -
\begin{bmatrix}
i_{\text{d load}} \\
i_{\text{q load}} \\
\end{bmatrix} + \text{C} \omega
\begin{bmatrix}
v_\text{q load} \\
-v_\text{d load} \\
\end{bmatrix}
$$
From there, I derive plant model and control structure in $dq$ reference frame (only $d$ frame is shown):

In this case disturbance  $i_{\text{d load}}$ is fed directly to the plant model. i.e. $P_3=1$. Open-loop inner loop equation is:
$$
G_\text{co}(s) = \underbrace{K_\text{P} \left( 1 + \frac{1}{s T_\text{I}} \right)} _\text{PI controller} \underbrace{ \left(\frac{1}{s \text{L}}\right) } _{\text{Process}}
$$
and closed loop ($P_1$) is
$$
G_\text{cc}(s) = \frac{G_\text{co}(s)}{1 + G_\text{co}(s)} = ... = \frac{s T_\text{I}+1}{s^2 \text{L} T_\text{I}/K_\text{P} + s T_\text{I} + 1}
$$
Compensator transfer fucntion is $G_{\text{ff}}=1/G_\text{cc}(s)$ and numerator is a higher order than the denominator.

I cannot implement this transfer function in MATLAB Simulink. Is there a way forward?
I don't see any issues with the method! Is this the right approach? Articles deal little with feed-forward signal or cross-coupling compensation.
When I make $P_3$ equal to $P_1$, i.e. $G_{\text{ff}}=1$, the disturbance is eliminated. In other cases, there is a huge overshoot for step disturbance. It eventually stabilises.
I also have to compensate the cross-coupling element $\omega\text{C}v_{\text{q load}}$ using the same compensator?
In the real system (e.g. DSP), I will implement only the control part in $dq$ frame. Do I need to add compensators to load current and cross-coupling signals? Also, the plant's transfer function will be different and unknown. How to determine compensator there or leave it at 1 and only tune PI controller to get the desired response.

1 - What transfer function are you trying to implement? The controller transfer function or the "plant" transfer function?
2 - You are trying to simulate some kind of inverter/active rectifier, right? Inverters are non-linear, you can simulate them with only linear blocks if you want but you will not be able to test your controller with the non-linearities that are present in real-life.
3 - What are P1 and P3 ?
4 - You kind of compensate your cross-coupling terms in the second image the $\omega Li_q$  and  $\omega Li_q$ terms. The cross-coupling terms are compensated to improve the response of your controller, but they are not technically needed.
5 - In a real-life implementation you don't have to compensate the cross-coupling terms but compensating the cross-coupling makes your control system settle faster if I recall correctly.


Polynomial division to separate FIR and IIR in cascade I have a FIR filter $L(z)$ and an IIR filter $H(z) = \frac{B(z)}{A(z)}$ which are cascaded together as $L(z)H(z) = L(z)\frac{B(z)}{A(z)}$.

For specific reasons I want to calculate a parallel topology that is equivalent to the original $L(z)H(z)$. I can achieve this using long division to divide $L(z)B(z)$ by $A(z)$ to obtain:

$L(z)H(z) = \frac{L(z)B(z)}{A(z)} = Q(z) + \frac{R(z)}{A(z)}$

Where $Q(z)$ is the quotient and $R(z)$ is the remainder of the long division of $L(z)B(z)$ by $A(z)$.

However, the resulting IIR, $\frac{R(z)}{A(z)}$, has a delay equal to the number of taps of the original FIR $L(z)$

Is it possible to do some clever re-arranging before the division so that the resulting IIR part of the long division is not delayed?

Of course one approach is to flip the filters so the magnitude of the poles are greater than one and perform long division using an unstable filter and flip back (which is a simple change of variable $z \to \frac{1}{z}$ and back again), but this isn't a useful solution.

The way that you currently split the FIR from the IIR part results in an FIR filter that exactly matches the first $N$ samples of the impulse response of the original cascade. Consequently, the parallel IIR filter must have $N$ initial zeros in its impulse response, otherwise it would mess with that first part of the original response which is already taken care of by the FIR filter. If you choose another FIR filter, then the corresponding IIR filter won't have as many initial zeros (i.e., delays without multipliers in its implementation). The question is why you would want this, because in that case the number of additions and multiplications necessary for implementing the IIR filter would increase. You would get the same number of delays, but instead of just delaying the input you would get non-zero filter coefficients.

Note that if $Q(z)$ denotes the transfer function of the FIR filter, and $R(z)$ is the numerator of the IIR filter, you just need to satisfy the following equation:

$$R(z)=L(z)B(z)-A(z)Q(z)\tag{1}$$

There are infinitely many solutions. You can choose any arbitrary FIR filter $Q(z)$, and the numerator $R(z)$ of the corresponding IIR filter follows from $(1)$. The most efficient solution, however, is obtained by choosing the FIR filter that is just a truncated version of the impulse response of the original cascade. In that case, the IIR filter $R(z)/A(z)$ delays the input by the length of the FIR filter, and takes care of the remaining tail of the original impulse response. That solution is just the one obtained by long division.

Also take a look at this related question and its answer.


What is the name for a constant-heading Kalman filter model for vehicle tracking? When applying Kalman filtering to estimate the position of a car, there are several different vehicle dynamics models that you could use.
One of the simplest is "constant velocity" or CV, which uses the state
$$
\mathbf x = [x, y, v_x, v_y]
$$
A more sophisticated model is "constant turn rate and velocity" or CTRV, which uses the state
$$
\mathbf x = [x, y, \theta, v, \omega]
$$
My coworker has implemented an algorithm that's somewhere in between these two, and I'm wondering what to call it. It uses the state
$$
\mathbf x = [x, y, \theta, v]
$$
I would expect that this model is called "constant heading and velocity" or "constant bearing and velocity" or similar but I haven't found any papers describing it.
Is there a name in the literature for this model?

I think the magic acronym is CHCV, "constant heading constant velocity". This returns at least a few results on Google.


Confusion in Peak to Average Power Ratio (PAPR) (1) What does it mean by 'a system with low peak to average power ratio requirement'? Doesn't every system idealises it to be lower? Or there are some exceptions to it?

(2) Higher PAPR can be a problem in signal carrier modulations too, I read the formula for single tone PAR, but I'm unable to visualise it.

Preference for low PAPR is subjective. From a power efficiency and unit cost perspective, low PAPR is certainly beneficial, but from a spectral efficiency perspective, it typically isn't. High entropy modulation means high PAPR (i.e. signal looks more like Gaussian random process).

A low PAPR can mean that less backoff is needed in a power amplifier, and hence better power efficiency can be achieved. Similarly, low PAPR means less headroom required in an ADC and thus allows for lower-precision in the ADC. These back-off/headroom requirements for TX and RX are closely related to PAPR. Some typical systems considerations that are linked to system requirements:

modulation - higher-order QAM and/or multicarrier modulation means higher PAPR, which means an increased budget for headroom across signal paths.
fading - not normally a part of PAPR, backoff for fading is closely related and arguably becomes an overall PAPR consideration for RX paths (we normally reference PAPR as signal-related, but fading typically increases the ratio)
code rate and blocksize - high coderate and large blocksize, particularly at high-order modulation, means that the signal can be highly susceptible to clipping, non-linearity, and other impairments, and thus we need to allow headroom for the rarest of peaks (definition of PAPR is somewhat subjective, too; are we worried about the 1E-3 probability peak or the 1E-6 probability peak?)
bitwidth - ADC, DAC, and processing chain bitwidth. Higher PAPR plays here, too. For higher PAPR, more backoff from full scale is needed to avoid distortion of the peaks. This requires more bits and analog linearity for signal representation and processing.

For multicarrier modulation (OFDM, DMT, etc), PAPR can certainly be problematic. To visualize, try modulating a high-order QAM constellation on subcarriers using an FFT (IFFT). The sum of tones results in a constructive and destructive combining across the IFFT symbol, giving high PAPR.


PSD and $\lim_{T\rightarrow \infty} \frac 1 {2T} \int_{-T}^T x(t)\bar y(t)\,dt$ From Wikipedia, I taken a definition of power spectral density:

For continued signals that describe, for example, stationary physical processes, it makes more sense to define a power spectral density (PSD), which describes how the power of a signal or time series is distributed over the different frequencies, as in the simple example given previously. Here, power can be the actual physical power, or more often, for convenience with abstract signals, can be defined as the squared value of the signal. For example, statisticians study the variance of a set of data, but because of the analogy with electrical signals, it is customary to refer to it as the power spectrum even when it is not, physically speaking, power. The average power P of a signal x(t) is the following time average:
$$P = \lim_{T\rightarrow \infty} \frac 1 {2T} \int_{-T}^T x(t)^2\,dt.$$

I premise that I'm not an expert on the signal theory, therefore I do apologise if this question isn't much precise.
Let $x(t),y(t)$ be two complex signals and denote with $\bar y$ the complex conjugate of $y$. What represents the following formula?
$$\lim_{T\rightarrow \infty} \frac 1 {2T} \int_{-T}^T x(t)\bar y(t)\,dt$$
Is it related to PSD?
Thanks in advance.

For power signals $x(t)$ and $y(t)$, the function

$$R_{xy}(\tau)=\lim_{T\rightarrow\infty}\frac{1}{2T}\int_{-T}^{T}x(t)\bar{y}(t+\tau)dt\tag{1}$$

is the cross-correlation of $x(t)$ and $y(t)$. So the expression you're asking about is the cross-correlation of $x(t)$ and $y(t)$ evaluated at lag $\tau=0$:

$$R_{xy}(0)=\lim_{T\rightarrow\infty}\frac{1}{2T}\int_{-T}^{T}x(t)\bar{y}(t)dt\tag{2}$$

Since the Fourier transform of $(1)$ is the cross-power spectrum $S_{xy}(\omega)$ of $x(t)$ and $y(t)$, the expression in $(2)$ is related to $S_{xy}(\omega)$ by:

$$R_{xy}(0)=\frac{1}{2\pi}\int_{-\infty}^{\infty}S_{xy}(\omega)d\omega\tag{3}$$

So the expression in your question is proportional to the integral of the cross-power spectrum of the power signals $x(t)$ and $y(t)$.


Absolutely summable signals Given an absolutely summable signal $x[n]$, the $z$-transform $X^z(z)$ is rational with a pole at $z=0.5$.

Given the following the statements:

$x[n]$ has a finite support in the time domain.
$x[n]$ is a left sided signal.
$x[n]$ is a right sided signal.
$x[n]$ is a two sided signal.

Can anyone please tell which statements are true and explain why or why not?

Absolute summability implies that the Fourier transform converges uniformly, so this means that the unit circle is contained in the region of convergence (ROC).

Now, the pole at 0.5 is the limit of the ROC, so the ROC is from 0.5 outwards (in order to include the unit circle). This means that the signal is right-sided.

Finite support signals cannot have poles (other than $z=0$ or $z=\infty$).


Why is it possible to frequency sweep with just one sinusoid? One way to measure the frequency response of a filter is to input separate cosine/sine wave at the desired frequency, measure the steady state output, stack the inputs/outputs into complex arrays, and divide. A more efficient way is to use the DFT, and somehow we get away with just inputting one sinusoid.

Input $x_n = A \cos(\omega_0 n)$, with $\phi = 0$
Measure the steady-state output:
$$
  \begin{align*}
    y_n =  A | H(e^{j \omega_0}) | \cos\left (\omega_0 n + \angle H(e^{j\omega_0}) \right )
  \end{align*}
$$
Compute the DFTs $Y_k$ and $X_k$
Record $H_{\hat{k}} = \frac{Y_{\hat{k}}}{X_{\hat{k}}}$ at $\hat{k}$ that corresponds to $\omega_0$

What I'm trying to figure out is why this works, clearly there's phase information but I'm missing something in the derivation.

$$
\begin{align*}
  X_{\hat{k}} &= \sum_{n=0}^{N-1} x_n e^{-j \omega_0 n} \quad \text{where } \frac{2\pi}{N} \hat{k} = \omega_0 \\
    &= \sum_{n=0}^{N-1} A \cos(\omega_0 n) e^{-j \omega_0 n} \\
    &= \frac{A}{2} \sum_{n=0}^{N-1} \left (e^{j\omega_0 n} + e^{-j\omega_0 n} \right ) e^{-j\omega_0 n} \\
    &= \frac{A}{2} \sum_{n=0}^{N-1} 1 + e^{-j 2 \omega_0 n}
\end{align*}
$$

$$
\begin{align*}
  Y_{\hat{k}} &= \sum_{n=0}^{N-1} A \left | H(e^{j\omega_0}) \right | \cos\left (\omega_0 n + \angle H(e^{j\omega_0}) \right ) e^{-j \omega_0} \\
    &= \frac{A}{2} \left | H(e^{j\omega_0}) \right | \sum_{n=0}^{N-1} \left (e^{j(\omega_0 n + \angle H(e^{j\omega_0}))} + e^{-j(\omega_0 n + \angle H(e^{j \omega_0}))} \right ) e^{-j\omega_0 n} \\
    &= \frac{A}{2} \left | H(e^{j\omega_0}) \right | \sum_{n=0}^{N-1} e^{j \angle H(e^{j\omega_0})} + e^{-j \angle H(e^{j \omega_0})} e^{-j 2 \omega_0 n} \\
    & \; \Big\Downarrow \; \text{?} \\
    &= \frac{A}{2} \left | H(e^{j\omega_0}) \right | e^{j \angle H(e^{j\omega_0})} \sum_{n=0}^{N-1} 1 + e^{-j 2 \omega_0 n}
\end{align*}
$$

The idea is that the phase information comes out by dividing:

$$
\begin{align*}
  H_{\hat{k}} &= \frac{Y_{\hat{k}}}{X_{\hat{k}}} = \left | H(e^{j\omega_0}) \right | e^{j \angle H(e^{j\omega_0})}
\end{align*}
$$

What am I missing?

Thanks to geometrikal and Jason R, the sum of exponentials is equal to zero when $\omega_0 = \frac{2\pi}{N} \hat{k}$ and $0 \leq \hat{k} \leq N-1$, which is obvious when you realize that it's just the sum over $2 \hat{k}$ periods of $e^{-j2\omega_0n}$, where $T = \frac{2\pi}{2\omega_0} = \frac{N \pi}{2\pi \hat{k}} = \frac{N}{2 \hat{k}}$.

For $\hat{k} > N-1$, I suppose you would have to zero pad, which is where the digital sinc Jason R. mentions comes in.


float image to integer type image I have a multispectral satellite image and I would like to "transform" each of the image bands into one with integer values. For example to stretch? the image between 0 and 255. Not sure how to proceed with this. I also thought I could make a histogram and replace each pixel value with the bin edge in which it falls into. I'm not sure. I'm currently working in python. Any thoughts on this?

You should check the Envi documentation (a standard tool for working with multi/hyper spectral images) on the stretches it uses to fit multi/hyperspectral data into 8 bits for display.

Here is the page on some of the stretches used, which apply an affine function to the pixel values and clip the values to the lowest and highest displayable value. In particular, look at the optimized linear stretch, which is the default one Envi uses to display things.


consecutive filtered blocks with filtfilt I'm doing basic signal processing in Python. I'm filtering consecutive blocks of a signal with scipy's filtfilt (to get zero phase). I would like to concatenate several blocks after filtering, however, each filtered block doesn't seem to be properly continuous with the other filtered blocks. Is there some way to work this out?

I'd like to note that when I'm doing the same with scipy's lfilter with "Initial Conditions"(zi) parameter being passed to it at each iteration the filtered blocks seems to be continuous.

Any help will be highly appreciated.

Mike

Both the SciPy and MATLAB implementations of filtfilt use a mirrored part of the beginning and end of the signal (roughly with a length equal to the impulse response of the filter) to reduce the effect of initial conditions on the filtered signal (i.e. to reduce ringing the signal is artificially lengthened by concatenating mirrored parts of the signal).

As these mirrored signal parts are probably not equal to the data in the previous and subsequent data blocks, concatenation of the filtered blocks will not result in a continuous signal.

A possible workaround might be to select longer data blocks from the original signal and discard a part from the beginning and end after application of the filter but before concatenation. Note that this might require a FIR filter to work properly although I would expect acceptable results for a IIR filter if the data block can be made long enough.


Understanding kurtogram parameters I am about to understand a kurtogram, and don't understand what means the value of "K" (presented in table 1),or especially why takes values of 1.6 ; 2.6 ; 3.6 etc.
Other question is how do interpret the kurtogram (figure 2)

[Given figures are from A new improved Kurtogram and its application to planetary gearbox degradation feature analysis, Xianglong Ni et al.]
This $k$ in the graph (and the $K$  in the Table) corresponds to an equivalent "dyadic-power", providing the number of subsequences at each level, alternating dyadic splits with filters on subbands $[0,1/4]/[1/4,1/2]$ and $[0,1/6]/[1/6,1/3]/[1/3,1/2]$.   The number of such subsequences in the 1/3-binary tree are $2,3,4,6,8,12\ldots$, approximately $2^1,2^{1.6},2^{2},2^{2.6},2^3,2^{3.6}\ldots$, since $\log(2)/\log(3)\sim 0.6$.
Additional sources:

Calcul Rapide du Kurtogramme et Applications, Jérôme Antoni (short, in french)
Fast computation of the kurtogram for the detection of transient faults, Mechanical Systems and Signal Processing, 2007

I unfortunately cannot help totally on the interpretation of this particular paper. I remember that the original  spectral kurtogram was meant to detect, localize or characterize nonstationarities from signals. The kurtogram was a fast version, reusing principles from wavelet packets and multiband filter banks and quasi-analytic filters.  When looking at maximum values in a kurtogram, if they are concentrated, one may suspect non-stationary components in a noisy signal. Then, it is possible to extract or invert the aforementioned component to check or denoise the original data.
I would suggest you to read several papers, especially the early ones, showing interpretations of different kurtograms.


What does the variable $\tau$ mean in frequency modulation? Frequency modulation formula:
$$
x_{FM}(t) = A_{c}\cos\left(2\pi f_{c} t+ 2\pi k_{f}\displaystyle\int_0^{t}m(\tau)d\tau+\phi_{0}\right)
$$
Its change from $m(t)$ to $m(\tau)$.
What is $\tau$? I wonder what it means.

It's just a dummy integration variable. The time variable $t$ is the upper integration limit, so you can't use it as the integration variable. An expression like
$\int_0^tm(t)dt\qquad(???)\tag{1}$
doesn't make sense, whereas
$$\int_0^tm(x)dx\tag{2}$$
does.
Take a look at this question and its answers over at math.SE.


How to Apply Least Squares Estimation for Sparse Coefficient Estimation? The model is expressed as, $$y(n) = \sum_{i=0}^{p-1} r(i) x(n-i) + v(n) \tag{1}$$ where  $\mathbf{r} = [r_1,r_2,\ldots,r_p]^T$ is the sparse channel coefficients of length $p$, $\mathbf{x} = [x_1,x_2,...,x_{n-p}]^T$ is the one dimensional input and v is additive White Gaussian noise of zero mean. y = filter(.) command is used to model the above equation and thus creating an FIR filter or a moving average (MA) model. The order of the MA model is p=3.

So, $y = [y(1),y(2),....,y(100)]$ is a vector of 100 elements.  I am generating noise of variance 0.1 I want to estimate the sparse channel coefficients using LASSO. As there are p channel coefficients, I should get p estimates.

According to the equation of LASSO, ||rx - y||_2^2 + lambda * ||r||_1 I am estimating the sparse coefficients, r. As the true coefficient array contains p elements, I should get p estimated elements. I am not quite sure if this is the way to do. I have not found any example on LASSO being applied to univariate time series model such as ARMA. I don't know how to estimate the sparse coefficients using the appropriate algorithm and need help.

The first part of the Equation : $||rx - y||_2^2$ is  a least squares formulation which I can solve using Least Squares Approach. In order to implement LS, I have to arrange the input in terms of regressors. However, if the coefficients, $\mathbf{r}$ are sparse then I should use LASSO approach. I have tried using Matlab's LASSO function. For LASSO, I rearranged the input data $x$ in terms of regressors, but I don't know if this the correct approach.

I need help. Is there an approach to include the sparsity term in the LS?

Please find below the code for LASSO using Matlab function. As a toy example I am just assuming model order to be of lag 3 but I know that LASSO can be applied efficiently to a large model. I can test for larger order MA model having lag > 3.

% Code for LASSO estimation technique for
%MA system, L = 3 is the order,

%Generate input
 x = -5:.1:5;

r = [1    0.0   0.0];% L elements of the channel coefficients
%Data preparation into regressors
X1 = [ ones(length(x),1) x' x']; %first column treated as all ones since    x_1=1

y = filter(r,1,x); % Generate the MA model
[r_hat_lasso, FitInfo] = lasso(X1, y, 'alpha', 1, 'Lambda', 1, 'Standardize', 1);

OUTPUT :

The estimates returned are r_hat_lasso = 0, 0.657002829714982, 0

Question : This differs very much from the actual r.
Is this wrong?

I am not entirely sure what matlab's LASSO routine does so I started with Ordinary Least Squares (OLS) and worked backwards.  From an OLS perspective X1 as you have it won't work.  You've got a regressor that is as all ones, but your parameter inputs for you example data (r) doesn't contain an offset.  Essentially your model doesn't fit you're data well, which is why you get a weird answer.

Second, you have repeated regressors which is going to confuse things even more.  I've rewritten your code to shift your regressors to line up the three inputs that would result from MA using p = 3.

You can perform regular old OLS on the data X1 and y to get the exact input r.  If I use matlab's LASSO function I can get the proper components of r but they are off by a scalar.

    % Code for LASSO estimation technique for
    %MA system, L = 3 is the order,

    %Generate input
    x = -5:.1:5;
    xm = circshift(x,[0,1]);
    xm2 = circshift(x,[0,2]);

    r = 0.5*[1.0    1.0   0.0]'; % L elements of the channel coefficients

    %Data preparation into regressors
    X1 = [ x' xm' xm2'];

    y = X1*r;
    y1 = filter(r,1,x); % Generate the MA model, so I can plot and check.

    r_hat = X1'*X1\(X1'*y);

    [r_hat_lasso, FitInfo] = lasso(X1, y, 'alpha', 1, 'Lambda', 1, 'Standardize', 1);

Output:

    r_hat =
    0.5000
    0.5000
    0.0000

    r_hat_lasso =
    0.3235
    0.3231
         0


For discrete wavelet transforms, is redundancy needed to ensure perfect reconstruction? So I was using pywt and I noticed that in the padding mode documentation they explain the following:

  DWT performed for these extension modes is slightly redundant, but ensures perfect reconstruction.

So I was wondering whether for discrete wavelet transform this redundancy was needed theoretically or if it was an implementation choice. This person at pywt seems to think it is theoretical but I wanted to know why.

Theoretically, the DWT (or the discrete wavelet transform) is (under some conditions) a discrete sampling of the continuous  wavelet transform allowing perfect synthesis, for "infinite length" discrete data. But this yields "discrete" infinitely many coefficients.
In some cases, There exist some wavelets, like Haar's, that can be implemented in a non expansive manner.

Practically, for discrete length signals that are not "even" (like with a power-of-two length), one must find options related to signal and wavelet symmetries. This is an important issue in data compression, as one often does not want to increase data size before entropy coding (but this is another debate).

To remain concrete, think about signal $[1,-1,3]$. How would you represent it with only 3 coefficients, knowing that the dyadic version should have the same number of approximation and detail coefficients? Many tricks have been developed for that. However, if you are not concerned with compression, most can afford a little redundancy.


How to properly match bit growth of a FIR filter with a reference value I am trying to implement a chain of CIC/FIR filters on an ZYNQ FPGA.
Using the Xilinx FIR compiler works fine so far but I am unable to properly get all the math.

At the moment I have 2 chains of filters:

$5\times$ downsampling, passband $12.5\text{MHz}$:
→ FIR filter 2, 5x downsampling, input 17.7.0, output 17.7
$625\times$ downsampling, passband $100\text{kHz}$:

→ CIC filter, $25\times$ downsampling, input 14 bit, output 16 bit

→ Compensation FIR filter, no downsampling, input 17.7, output 17.7

→ FIR filter 1, $5\times$ downsampling, input 17.7.0, output 17.7

→ FIR filter 2, $5\times$ downsampling, input 17.7.0, output 17.7

I have a 14 bit, $20\text{kHz}$ signal that I feed into both chains:

If I feed it into the first chain I sign extend the 14bit value to 17bits and add 7 zero-bits at the end (.7).
I get a nice output that wins me one bit.
If I feed the signal into the second chain I directly feed the 14bit signal into the CIC filter and left-shift the output signal to only keep the 16 MSBs. I then sign extend that 16bit value to 17bit (yes I could have also just kept the first 17bits of the CIC return value, if that would be better, please tell me) and add 7 zero-bits at the end (.7) as well.
The filters win me 1-2 bit until the output of the "FIR filter 1" but at the output of the "FIR filter 2" I have an output where I lost 8bit (upper 8 non-sign bits are 0). This is the exact same FIR filter that gives me +1 bit in the first ($5\times$ downsampling) chain ... This broken chain also seems to have a lot of jitter after one period ... (I can exclude overflows as the source of error as the Xilinx FIR compiler guarantees no overflows).

How can this be, what am I missing?

In the same question I would like to ask:

Let's assume all the chains finally work properly. I then will have different bit growth. Initially I knew that 14bit correspond to $2.2\text{Vpp}$ which means that 12.86bits correspond to $1\text{V}$. But now I suddenly have a 17.7bit output. How do I know where my $1\text{V}$ reference level is now? Because it has shifted some.
I am quite sure that accounting just for the DC gain yields wrong results.

My waveforms:

Output $5\times$ chain:

Chain $625\times$:

I am very grateful for any hints and will gladly add more information if some more is required.

I've worked a lot with FIR and IIR filters in FPGAs, here are a few pointers

1 - Keep track of your bits and their respective weights.

2 - Keep track of your filter gain and the number of bits for the coefficients. For example, assume you have a 14-bit signed input. Your input can go from +16383 to -16384. Now you filter it in a 17-taps FIR low-pass filter (for example), its DC gain should be '1'. Assuming 16-bit coefficients, and no coefficient greater than or equal to '1', the sum of the coefficients should be 32768 (2^15). It could be slightly smaller, it doesn't really matter, especially if you round the output. However, it is possible, that because of the scaling, the sum of the coefficients could be 32769, 32770, etc. In this case your DC gain would be slightly bigger than unity and you will likely need guard bits.

3 - Assuming there can be no overflow at the output of the filters, now you need to determine how many "fractional bits" you need to keep. There are no simple answers. Assuming no decimation, it could be a good idea to keep the output format the same as the input. For example, 14-bit input with 1 sign-bit ---> 14-bit ouput with 1 sign bit. However, please round the data before truncating, otherwise you could insert some DC offset or other unpleasantness.

4 - If you decimate, one of the reasons why could be to increase the resolution of your data, in this case your output should have more fractional bits than your input.

5 - If your filters have a gain higher than 1 (on purpose, not because of quantization), such a resonant IIR filter, you will need to implement the appropriate number of guard bits. For example, you add a 20 dB-gain at a certain frequency, you will need 4 additional guard bits at the output.


Non-polynomial Z-transform My prof said that when a transfer function described by a z-transform is not polynomial, then i can't perform the anti-transformation.
But, what does it means to be not polynomial ?
Can you explain to me ?

First of all - I myself am not a pro in Control theory, but a mathematician - so I write what I think might be what your prof means.

Part where i am quite certain:

Not to be polynomial means, that there exists NO polynomial which describes the transfer function. E.g. let $H$ be your transfer function, then there exists no $n\in\mathbb{N}$ and $p\in P^n[X]$ such that $$H(s) = \sum_{i = 0}^{n}p_is^i$$
holds.
Examples for functions which are not polynomial are: exponential functions (e.g $e^s$) or trigonometric functions (e.g. $sin(s)$), but there are many more.

The next part is where I am not 100% sure, because it is not my main subject.

I assume you mean the inverse Z-Transform when saying anti transform. I do not know too much about the z-transform, but i think it is the discrete-time equivalent to the laplace-transform and therefore:

Let $H$ be a transfer function. Then $H$ is the Laplace-transform of an impulse response. Thus the inverse Laplace transform (i.e. $h$) exists. (The stability however is something different.)

But if you want to check for a given "transfer function" (an arbitrary function in the "s"-Domain) if there exists an impulse response corresponding to this "transfer function", then I think you need the inverse Laplace transform to converge, i.e.
$$ \frac{1}{2\pi i} \lim_{y\rightarrow\infty}\int_{\gamma-iy}^{\gamma+iy}e^{st}H(s)\,ds$$ to exist for some $\gamma\in\mathbb{R}$.
The integral then is your impulse answer. There exist non-polynomial functions for which this integral exists. Thus, there are non-polynomial transfer functions which correspond to a system. Many many authors in literature however assume only rational valued transfer functions... i really do not know why.

If this answer is NOT correct (because for the discrete-time case there is something mayor different), or does not help, please say so and apology. As I said, this is not my main subject.


Gradient Domain Reconstruction - Scaling Problem I am implementing reconstruction of image from gradient domain. This requires solving the following partial differential equation (a Poisson equation) on a 2D grid:

$$\nabla^{2}I=\mathbb{div} G$$

$\mathbb{div} G$ is divergence of a gradient field, which is simply an image with Laplacian operator applied (second-order derivative):

The original image have range of values 0 to 255 (8-bit image). The divergence image has values ranging from -341 to +318, which is expected.

However, when I solve the above equation for $I$, the result have very small scale. The values of reconstructed image range from -0.000038 to +0.000029.

I know the gradient domain stores information up to a constant. In this case however, the scale has been changed considerably.

When I stretch the reconstructed image to 0-255 interval, the correct result is obtained:

I am using Poisson solver based on the Full Multigrid Method (FMG) which works spectacularly in different application (HDR images). Here the values are stretched anyway so I never noticed the huge change of scale.

Is this scale change normal for gradient domain reconstruction or may there be

My ideas:

Some glitch in the numerical solution.
The $\mathbb{div} G$ image is padded with zeros so its size conforms with the solver - could this affect the scale?

The Poisson equation may suffer from ill-posedness if the boundary condition is not suﬃcient to yield a unique solution to the problem. The most often used boundary condition is Dirichlet boundary condition, in which the values of $I(x,y)$ are specified for those points $(x,y)$ on the boundary of $I$.

In your case, after you stretch your solution to another scale, you got the correct solution. It is possible the solution to your Poisson equation is not unique, and your algorithm just found out the correct one based on some restrained part (such as minimal norm of least square) in the optimization. It is much much more possible that it is already sufficient to determine $I(x,y)$ up to an arbitrary additive constant, and this constant is determined by the solution of Laplacian equation: $$\nabla^{2}\hat I= 0$$


Signal Reconstruction in Compressed Sensing with a Simple Vector Signal as an Example While going through the different types of reconstruction algorithm as mentioned in Richard G. Baraniuk - Compressive Sensing - Lecture Notes (Also on DocDroid), I came to know that minimum $ {L}_{1} $ norm reconstruction is better than minimum $ {L}_{0} $ and $ {L}_{2} $ norm reconstruction.

Can anyone show me why $ {L}_{1} $ norm is better by taking a simple signal vector of small length as an example?
I would like to know what will be the final optimization results in each reconstruction method.

Better is hard to measure unless there is a well defined measure.

Yet in the classic Compressed Sensing model (Sparse Representation) the Holly Grail is actually the $ {L}_{0} $ Pseudo Norm.

The problem is $ {L}_{0} $ isn't really a norm and hence the optimization problem isn't convex.
Moreover, it is not smooth (It has discrete properties) and the problem:

$$ \arg \min_{x} \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2} \quad \text{s.t.} \quad {\left\| x \right\|}_{0} \leq k $$

Is actually NP hard. Namely we don't know how to solve it efficiently.

We can prove that under some cases solving the following problem:

$$ \arg \min_{x} \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2} \quad \text{s.t.} \quad {\left\| x \right\|}_{1} \leq k $$

Yields the same solution. Yet this form is Convex and we have many effective tools to solve it.

Hence it is preferable to use this model in the cases it is good enough.


Meaning of concentric circles on the spectrogram I've got this unusual spectrogram and I am really curious what these circles could mean. Has anyone ever seen anything like this? What could be the source of the signal?

In Understanding FFT Overlap Processing, p. 10, it is suggested that uniform of constant windows, for certain hop patterns, with repeated pulses, could generate such artifacts.


Modulation and allocation in GSM900 I want to do one project and I just want to make sure that I've understood correctly.

Human's voice frequency is around 85-255Hz.

When we talk on phone, carrier is modulated with voice signal in order to shift frequency spectrum to higher frequency for transmission. If we take for example GSM900 which is 876 - 915 MHz paired with 921 - 960 MHz. That means that our speech will be in this GSM900 range after modulation?

Moreover, how much bandwidth user gets in 876 - 915 MHz paired with 921 - 960 MHz range? Any articles about that?

GSM since it's first day has been digital (Fat32 is incorrect). That means you can't send analog audio. Luckily, analog audio had already been phased out in the 80s, in telecom.

In digital telephony, audio was digitally sampled at 8 kHz, 8 bits per sample, for a 64 kbit/s signal. If you has an ISDN phone, you got this quality at home, but it certainly was used inside and between the central phone switches.

GSM built on that 64 kbit/s idea, but used audio compression. The first codec used was GSM Full Rate, which achieved 5x compression for a 13 kbit/s net data stream.

GSM furthermore divides every frequency band into 8 slots, which means that phones can only transmit short 5 ms bursts and then have to wait for 7 other bursts from other phone users. Band spacing is a mere 200 kHz; a typical carrier gets 5Mhz blocks of spectrum which means 25 bands and 200 slots total.


I somehow "proved" that given any LTI system, its transfer function has to be constant. What am I missing? A transfer function is defined as the Laplace transform of the ratio of output to input.

Also, every LTI system has an eigenfunction. Given such eigenfunction as an input, the ratio of the output to input is a constant, which means the transfer function is constant.

When the input to an LTI system is an eigenfunction, then:

$$ x(t) = e^{j\omega_0 t} => y(t) = H(\omega_0) e^{j\omega_0 t}$$

where $H(\omega_0)$ is the associated eigenvalue of the LTI system; the transfer function evaluated at the particular frequency:

$$ H(\omega_0) = H(s)|_{s = j\omega_0}. $$

Here not the transfer function $H(s)$, but just one value of it is found by just one eigenfunction.

But if you use the general expression for the eigenfunction, then the input is $x(t) = e^{j\omega t}$ and the output will be $$y(t) = H(j\omega) e^{j \omega t}. $$ Here you assume that $H(j\omega) = a$ ; a constant,  and therefore in the fraction expression you will find out that $H(s) = \frac{Y(s)}{X(s)} = a = H(j\omega)$. Here, the constant $a$ is actually the unknown variable $H(j\omega)$...

Note that the eigenvalue $H(\omega)$ will be a complex-constant for a given particular frequency $\omega_0$. But as the given frequency changes, the constant (in general) will change too. It's just like an ordinary function $f(x)$ treated as a single value for a particular $x = x_0$ but a rather a varying mapping for general $x$ values...


Very fundamental question: by what number do you divide your ADC counts? What number is it correct to divide by your ADC counts as a first step in converting counts to voltage? For example, in a 10 bit ADC, is it correct to divide by 1024 or 1023?

I have always divided by 1024, which means that a full count would yield less than full scale voltage reading. I think this is a more accurate representation of the hardware, though, since there is usually an LSB of error. This is also what I've read in every datasheet, as well as various DSP literature and white papers from ADI and TI. Please prove me right (or wrong) with the detailed reasoning why!

Links to lit:
http://www.ti.com/general/docs/litabsmultiplefilelist.tsp?literatureNumber=slyw038b

http://www.analog.com/en/education/education-library/data-conversion-handbook.html

http://www.analog.com/en/education/education-library/scientist_engineers_guide.html

The correct answer depends on implementation and what you decide to call Full Scale; if the end points are included in the conversion mapping used, then the answer would be $2^n-1$ as I depict in the graphic below.

The easy way to see this is to consider a simple 3 bit converter which has 8 levels, for the conversion of a unipolar signal from 0 to 1V, where we specifically have mapped (in the design of the converter) 0V (+/- 1/2 lsb) to 000 and 1V (+/- 1/2 lsb) to 111:

Here we see clearly how with 8 levels the resulting mapping is given as

$$V_{out} = 1V\frac{b}{2^3-1}$$

where b is the digitized value of interest, in this case ranging from 0 to 7.

Because our digital range is inclusive of the end points of our voltage range, it would always be 1 less than the number of levels in this case.

Although not an ADC, but a counter example I could find showing the implementation dependence is A DAC implemented with an R-2R ladder (https://en.wikipedia.org/wiki/Resistor_ladder). This shows a mapping where a division of $2^n$ is used, and the full scale voltage IF based on the reference voltage used is never actually reached as an output of the DAC (but why not simply call "Full Scale" the voltage corresponding to the maximum digital value? - Here you see that it is simply a matter of how you decide to scale it).


Avoiding DC spike while decoding FSK signal I am using the Quadrature Demod module in GNU Radio to decode an FSK signal, but the signal is appearing quite noisy. It has two distinct peaks centered around the advertised frequency, but they are close together and only a little above the DC spike at the center.

In my first attempt, I set the baseband frequency of my HackRF SDR to the exact center frequency of transmission with a sample rate of $2\textrm{ MHz}$. I then applied a $30\textrm{ kHz}$ low-pass filter to the signal before doing the quadrature demod on the signal. What I got out was a near quiet signal during idle periods and something resembling digital data during the transmission, but a bit on the noisy side. It was easy to set a trigger and separate out individual packets, but hard to identify the bits. A screen shot of the start of a packet is here:

So I decided to try capturing offset from the center frequency by about $900\textrm{ kHz}$ and multiply the data by a corresponding complex signal to shift it back down to baseband. The signal now is well above the noise floor during transmission, however, my quadrature demod is capturing large amounts of noise during the idle periods. If I stop the capture at just the right moment, I can see what looks like a very clean data signal, but considerably weaker amplitude than the surrounding noise.

How do I get the best of both worlds so I can cleanly detect packets and decode the bits with GNU Radio?

My flow graph from the last test:

That looks very nice:

I don't 100% agree with your phrasing on this one:

  however, my quadrature demod is capturing large amounts of noise during the idle periods.

Well, the quadrature demod takes successive samples and gives you the phase advance between them. If you only capture silence, or more precisely, the noise between packets, the phase of each sample should (in theory, at least) be absolutely uniformly random – and hence, your phase increases, too, should be random. Remember: the amplitude of an FM-demodulated (analog or digital) signal is not proportional to the original signal's power at the receiver (that's the beauty of FM).

There's different ways to tackle this; a simple one is in the "Power Squelch" blocks in GNU Radio (which really do what their name suggests – make sure the signal is suppressed as long as its samples are below a certain |·|² threshold). You could add one before your quadrature demod.

Another method: The mean of your "demodulated noise" should be 0 – so if you low-pass filter your demodulated signal so that you can still make out the symbols, but the phase-noise induced "jumping" gets averaged to a 0, you could just have a decider afterwards that says "ok, this is a 1 or a -1 in a packet, or this is a 0 belonging to inter-packet time".

 Remarks:

Instead of multiplying with a complex cosine, you can use the "Rotator" block; should be a bit faster (not that it matters at these rates at all).
Even more comfortably: You can simply replace your low-pass filter, multiplier and cosine by a single "Frequency Xlating FIR". Use the "Low Pass Filter Taps" block to calculate the taps for the filter you want, and use the id of that block in the xlating FIR's "taps" field.
As noted in the comments: I recommend new designs use Qt instead of WX GUI, because WX is no longer maintained and will go away in future releases (i.e. GNU Radio 3.8)
The hard part might be recovering the timing of the symbols – Look at the synchronizers GNU Radio offers. Classically, many designs used Müller&Muller, but to be honest, that is an algorithm designed for implementation in hardware, not software. It's a bit involved to parameterize that, the "PFB clock sync" should behave better (it does, in fact, work similarly to M&M, but is quite a bit more accurate in approximating the derivative of the signal)
A quadrature demod is not the only thing that can decode an FSK: try filters that select each FSK frequency, followed by a complex to mag square each, and a "max" block. If you're clever about designing your filters (ie. don't overspecify things), this is, on PC-style hardware, not necessarily harder on the CPU than the $\arctan$ that quadrature demod has to do.


Amplitude after Fourier transform How to obtain the correct amplitude after the numerical Fourier transform of a signal?
Example: consider an exponential decaying wave $y(x)=e^{-x}\sin(100\pi x)$ with Fourier transform $y_f(x_f)$ which should be a peak with some broadening. I am searching for the appropriate approach to numerically approximate the continuous function $y_f(x_f)$ as close as possible. The calculation of the approximation of the function $y_f(x_f)$ should thus be independent of the number of samples and the range that y(x_i) is defined.
However, when I calculated the Fourier transform for several x-ranges (N in the code), I obtain different results:
fig, ax = plt.subplots()
N=[600,1200,2400,4800]
for i in N:
    N = i
    # sample spacing
    T = 1.0 / 800.0
    x = np.linspace(0.0, N*T, N)
    y = np.exp(-x/1)*np.sin(50 * 2.0*np.pi*x)

    #Fourier
    yf = scipy.fftpack.fft(y)
    xf = np.linspace(0.0, 1.0/(2.0*T), N/2)
    yp = 2.0/N * np.abs(yf[:N//2])

    label = "N="+str(N)
    ax.plot(xf, yp, label=label)

ax.legend()
plt.show()

The result gives different peak amplitudes for the different N-cases. The more sample points are taken (higher N), the lower the fft peak value... So, if you take N really high, then the Fourier transform amplitude decreases?
I would expect that the peaks have the same amplitude because $y(x)$ remains the same function in all cases. If the function does not change, the Fourier transform does also not change...
How to obtain the correct peak amplitude, such that the peak value is the same in the three cases?

Relating the DFT to the FT (CTFT) is a big issue.  Let's start with the basic definitions without any domain specification using a $\frac{1}{N}$ normalization on the DFT.
$$ FT(x(t))(f) = \int x(t) e^{-i2\pi t f } dt  $$
$$ DFT(x[n])(k) = \frac{1}{N} \sum x[n] e^{-i 2\pi \frac{n}{N} k } $$
They are obviously similar.  Assuming they cover the same interval, we can map the domain variables to each other.
$$
\begin{aligned}
t &= \frac{n}{N} = \frac{1}{N} \cdot n \\
f &= k  \\
( \Delta t  &=\frac{1}{N} ) \to dt
\end{aligned}
$$
[Edit: Does anybody want to argue the (1/N) doesn't belong in the DFT definition?]
The units of $f$ and $k$ are cycles per the interval.
Of course, we are assuming the sequence of values matches the function at every point.
$$ x[n] = x(t) $$
The FT is defined on a domian of $-\infty$ to $\infty$, and the the DFT is defined on a domain of $0$ to $N-1$.  The corresponding $t$ values are $0$ and $\frac{N-1}{N}$.  We can now include the domains in the definitions for a direct comparison on the same interval.
$$ FT_{DFT}(x(t))(f) = \int_{0}^{\frac{N-1}{N}} x(t) e^{-i2\pi t f } dt  $$
$$ DFT(x[n])(k) = \frac{1}{N} \sum_{n=0}^{N-1} x[n] e^{-i 2\pi \frac{n}{N} k } $$
To make the $FT_{DFT}$ transform fit the $FT$ transform, we have to introduce a rectangular window function on the interval.
$$ w(t) =
\begin{cases}
= 0 & t < 0 \\
= 1 & 0 \le t < 1 \\
= 0 & 1 \le t \\
\end{cases}
$$
$$ FT_{DFT}(x(t))(f) = \int_{-\infty}^{\infty} w(t) x(t) e^{-i2\pi t f } dt  $$
Please notice that there is no window function in the DFT! The window function is in the FT.
From these definitions, it is very clear that the FT is the limit of the DFT with a $\frac{1}{N}$ normalization factor as $N$ goes to infinity.
$$
\begin{aligned}
\lim_{N \to \infty} DFT(x[n]) &= FT_{DFT}(x(t)) \\
\lim_{N \to \infty} \frac{1}{N} \sum_{n=0}^{N-1} x[n] e^{-i 2\pi \frac{n}{N} k } &= \int_{-\infty}^{\infty} w(t) x(t) e^{-i2\pi t f } dt \\
\end{aligned}
$$
This is the "right" way to view it.  An integral is the limit of a summation.
It is possible to to view the DFT as a sampled version of the FT by using a second window "function" (loosely speaking) of a train of Dirac deltas.  But this is not the definition of the DFT, requires advanced Real Analysis understanding, and you are merely getting back to the definition by working backwards.
Some EE professors neglect to mention this, some just don't know it.
So, just like you can approximate an intergral with a summation, you can approximate an FT on a finite interval with a $\frac{1}{N}$ normalized DFT.
To get intervals other than $[0,1)$ on $t$, call your new time $\tau$ and rescale it.
$$ t = \frac{\tau}{T} = \frac{n}{N} $$
Feel free to shift it too.
$$ t = \frac{\tau- \tau_0}{T} = \frac{n}{N} $$
This makes the sequence definition equivalent to
$$ x[n] = x_{\tau}(\tau) = x_{\tau}\left( \tau_0 + \frac{T}{N}n \right) $$
The DFT does not care what your sampling rate is. aka $\frac{N}{T}$ with these names nor where you place the frame, $\tau_0$.

Followup:
Questions from the OP's comment under the post question:
(All answers assume you are dealing with a pure real tone with a slighlty varying envelope.)
"How do you get rid of the influence of the number of samples and the range over which y(x) is defined? "
By framing your interval on a whole number of cycles.  This is the principal justification in my mind for a $\frac{1}{N}$ normalization factor of the DFT.  On a whole number of cycles (for a pure steady tone) a single conjugate pair of bins are non-zero.  The value in the bin is independent of N when normalized.  Thus the reading you get matches the continuous case exactly and is independent of the sample count.  Also, there is no need to calculate the whole DFT (take a FFT), a single bin calculation will do.  With a consistent frequency having a set of basis vectors is most efficient.  With a varying frequency, a Goertzel calculation yields the equivalent more efficiently.
"How do you obtain the correct amplitude of the Fourier transform of the function y(x)? "
With a single bin, the magnitude of the bin is half the amplitude of the tone.  A small DFT frame size of 2 to 3 cycles in length ensures that the reading is fairly accurate on that interval.  The reading will be that of the best fit sinusoid on that frame.
"How do you properly take into account the number of samples and the range over which y(x) is defined? "
Not sure what you mean by this.  The sample rate and signal frequency define how many samples per cycle you will have.  Decide the number of cycles you want, multiply that by the number of samples in a cycle and you get the sample count (N) for that frame.  If you are calculating a single bin, there is no special advantage to powers of two or any other factor.


Checking if a lecturer can be heard clearly from the back of a room I would like to indicate to a human speaker standing at the front of a classroom if he/she is speaking loudly enough to be heard from the back. From some research, i have gathered that I am looking to calculate the Signal to Noise Ratio or something similar. It looks like this sort of thing can be calculated in MATLAB or preferably Python.

What types of inputs are needed? Is it possible to figure this out all from a mic in the back of the room? Or do I need a reference signal. For example, if i have a mic in the back of the room (connected via Bluetooth or WiFi) getting a noisy signal, will I also need a mic on (or near) the speaker so that i have something to compare to?

In the end, this would be running on a raspberry pi. Can anyone point me towards helpful libraries, modules, or tools for this type of thing?

This seemed like the best forum to ask this question. If there is a better place for it, let me know.

There's possibly an easy way out: Add two mics, one near the speaker and one in the back of the room. Run both signals through speech recognition software. If there is textual output from both and they agree, then the voice must be loud and clear. The software must be able to handle the language spoken though.

Other than that, there is research and standards on speech intelligibility. Perhaps it is better to start with a simplified model of intelligibility. Two important factors are speech level and noise level. By carrying out preliminary listening tests and measurements in the space, you can determine what the baseline noise level ($N_0$ dB) in the back is and how loud the speech must be ($S_0$ dB) in the minimum at the front mic for it to be intelligible in the back. Then during a lecture, you would detect voiced segments of the audio with the front mic and measure the speech level ($S$ dB) at the front during voiced segments and noise level ($N$ dB) in the back during non-voiced segments. If $S-N \ge S_0-N_0$, the speech should be intelligible. When measuring levels, you should weight the spectrum (by running the signal through a filter) so that level best reflects speech intelligibility or the detrimental quality of the noise. You would need to look into voiced/non-voiced detection, or simply have a threshold level for $S$ above which you classify the input as voiced.


how to detect two parallel lines in a binary image with MATLAB I am working on license plate detection using hough transform. I am trying to detect the horizontal and vertical edges and then to find the parallel lines.

My code:

img = rgb2gray(img);

img = medfilt2(img);

v = edge(img,'sobel','vertical');
v = bwlabel(v);
stats = regionprops(v, 'Area','BoundingBox','Image');
ids = find([stats.Area] > 30 & [stats.Area] < 2000);
v = ismember(v,ids);

Output:

Now I want to detect the parallel lines. I have consulted different examples on Mathworks but unfortunately I couldn't get much of it. Someone please guide me how to detect two almost straight and parallel lines in a binary image?

You mentioned the Hough transform, but your code doesn't use it. However, it can help you to find out the orientations of lines.

The maximum values of the Hough transform correspond to probable lines. These maxima are defined by their coordinates $(\rho, \theta)$, where $\rho$ is the distance from the origin to the closest point on the straight line, and $\theta$ is the angle between the $x$ axis and the line connecting the origin with that closest point (see Wikipedia).

Therefore, two parallel lines have a close $\theta$ (or $\theta + \pi$).


FFT method input argument have to be $2^n$ ? Does FFT method input argument have to be power of 2, i.e, $2^n$
I just realized there are many algorithm for FFT implementation, is there any algorithm that takes arbitrary amount of samples as input argument?
if yes, what is their advantages amongst each other?
namely, what algorithm does Matlab apply?

It is a common misconception that all FFT algorithms require their inputs to have a length that is a power of 2. The radix-2 Cooley-Tukey algorithm, the first to be widely known, has this limitation, but there are mixed-radix versions of the algorithm that don't require a power-of-2-sized input. In addition to Cooley-Tukey, there are a number of other FFT algorithms that have differing requirements on their input size; there are even techniques like Rader's algorithm that work on inputs whose sizes are prime.

Good, general purpose FFT implementations like MATLAB's fft() or FFTW support a wide range of algorithms and select one that is most appropriate for the specific situation. For more details on the specifics, go to the MATLAB documentation itself:

> doc fft

  Algorithms

  The FFT functions (fft, fft2, fftn, ifft, ifft2, ifftn) are based on a
  library called FFTW [3],[4]. To compute an N-point DFT when N is
  composite (that is, when N = N1N2), the FFTW library decomposes the
  problem using the Cooley-Tukey algorithm [1], which first computes N1
  transforms of size N2, and then computes N2 transforms of size N1. The
  decomposition is applied recursively to both the N1- and N2-point DFTs
  until the problem can be solved using one of several machine-generated
  fixed-size "codelets." The codelets in turn use several algorithms in
  combination, including a variation of Cooley-Tukey [5], a prime factor
  algorithm [6], and a split-radix algorithm [2]. The particular
  factorization of N is chosen heuristically.

  When N is a prime number, the FFTW library first decomposes an N-point
  problem into three (N – 1)-point problems using Rader's algorithm [7].
  It then uses the Cooley-Tukey decomposition described above to compute
  the (N – 1)-point DFTs.

  For most N, real-input DFTs require roughly half the computation time
  of complex-input DFTs. However, when N has large prime factors, there
  is little or no speed difference.

  The execution time for fft depends on the length of the transform. It
  is fastest for powers of two. It is almost as fast for lengths that
  have only small prime factors. It is typically several times slower
  for lengths that are prime or which have large prime factors.

    Note   You might be able to increase the speed of fft using the
    utility function fftw, which controls the optimization of the
    algorithm used to compute an FFT of a particular size and dimension.

  References

  [1] Cooley, J. W. and J. W. Tukey, "An Algorithm for the Machine
  Computation of the Complex Fourier Series,"Mathematics of Computation,
  Vol. 19, April 1965, pp. 297-301.

  [2] Duhamel, P. and M. Vetterli, "Fast Fourier Transforms: A Tutorial
  Review and a State of the Art," Signal Processing, Vol. 19, April
  1990, pp. 259-299.

  [3] FFTW (http://www.fftw.org)

  [4] Frigo, M. and S. G. Johnson, "FFTW: An Adaptive Software
  Architecture for the FFT,"Proceedings of the International Conference
  on Acoustics, Speech, and Signal Processing, Vol. 3, 1998, pp.
  1381-1384.

  [5] Oppenheim, A. V. and R. W. Schafer, Discrete-Time Signal
  Processing, Prentice-Hall, 1989, p. 611.

  [6] Oppenheim, A. V. and R. W. Schafer, Discrete-Time Signal
  Processing, Prentice-Hall, 1989, p. 619.

  [7] Rader, C. M., "Discrete Fourier Transforms when the Number of Data
  Samples Is Prime," Proceedings of the IEEE, Vol. 56, June 1968, pp.
  1107-1108.


How is a signal multiplied with the window functions? How the input signal get multiplied with the window functions like Hamming, Hanning, Blackman, etc?

I have tried to multiply the input signal of $10\textrm{ MHz}$ sampled at $100\textrm{ MHz}$ using the Blackmann window function. The output graph is

How is the multiplication perfomed ?

The input signal is $1\times 64$ row matrix, the Blackman window generates its samples in $64\times 1$ column matrix. The multiplication generates a $64\times 64$ matrix i.e one sample of input signal is multiplied to 64 samples of the window function. For a particular point sample there are 64 sample value. I intend to calculate 128-point FFT with padding zeros.

So for a particular sample point which value should be considered to calculate FFT?

Here is my code :

N = 64; % signal length (power of 2)
T = 10*(10^-9) ; % sampling period (and rate) is set to 1
A = 1; % sinusoid amplitude
phi = 0; % phase of zero
f = 10*(10^6); % frequency (under Nyquist limit)
nT = [0:N-1]*T; % discrete time axis
w = blackman(N,'periodic');
x1 = A*cos(2*pi*f*nT);
x = w*x1  ; % windowed signal
subplot (311);
stem (x1);
wp = [w;zeros((zpf-1)*N,1)];
subplot (313);
stem (x);
subplot (312);
stem (w);
disp (x1);
disp (w);
disp (x);

Matlab window functions such as hamming(N) returns the $N$ point Hamming window sequence $w[n]$ in a $N \times 1$ matrix (a column vector).

Therefore if your sequence $x[n]$ is represented in a row vector of size $1 \times N$ then for obtaining the sample-by-sample product $x[n]w[n]$, you should either transpose $w[n]$ or the other $x[n]$, so that their dimensions permit the correct sample-by-sample multiplication in between and produce the output sequence $v[n]=w[n]x[n]$ represented in a matrix of the size $N \times 1$ or $1 \times N$ depending on which one you have transposed.

v0 = x.*w      % x: Nx1, w: Nx1, v0: Nx1
v1 = x .* w'   % x: 1xN, w: Nx1, transpose w => v1: 1xN
v3 = x' .* w   % x: 1xN, w: Nx1, transpose x => v3: Nx1

Note that if the transposed signal is complex valued, then you should use the conjugate transpose instead :

v1 = conj(x').* w  % x[n] was complex valued...

Note also that the sample by sample product is perfomed by the $.*$ operation under Matlab;

  v = x .* w    % is a sample by sample product implementing v[n] = x[n]w[n]

If you omit the period and use the $*$ alone, you get matrix product or the dot product in case of vectors involved.

For this case when you have two vectors, $x$ of size $N \times 1$ and $w$ of size $1 \times N$, then their matrix product would produce a matrix of size $N \times N$

v = x * w    % x: N x 1, w: 1 x N => v: N x N


Multiply signal $x[k]$ with $\cos(2\pi\nu_0k)$, then given $X(\nu)$ draw resulting function in frequency domain? Let
$$y[k]=x[k]\cdot \cos(2\pi\nu_0k) .\tag{1}$$

Then, given a signal $x[k]$ with the DTFT $X(\nu)$ according to the following figure

what will the frequency domain for $Y(\nu)$ look like for a given $\nu_o$ in the range $0<\nu_0<a/2$?

According to my calculation, the resulting function $Y(\nu)$ is

$$Y(\nu)= \frac{1}{2}\left ( X(\nu-\nu_0) +X(\nu+\nu_0)  \right ). \tag{2}$$

This since the DTFT of cosine is two unit impulses in $\pm\nu_0$.

Is this correct? And how can this resulting $Y(\nu)$ be drawn in the frequency domain? Thank you! Some explanation with mathematical elements would be appreciated.

Would the following figure for $Y(\nu)$ be correct?

Thinking in terms of convolution with shifted impulses help. Multiplication in the time domain corresponds to convolution in the frequency domain. Your example is a classic showcase of frequency/band shifting using a "carrier" such as a cosine (or sine) or a complex exponential. Now, note, that the Fourier transform of these carriers are actually shifted impulses. When you convolve a function x(t) with a shifted impulse, say, D(t-T), the result that pops out is x(t-T), i.e., basically a shifted form of x(t). Note that the amount of shift is exactly equal to the shift in the impulse (the actual math is quite easy). Likewise, in the frequency domain, convolution with a shifted impulse results in your Fourier transform of the original signal getting shifted accordingly - For an impulse shifted to the right of origin, it feels as if the impulse "pushes all the bands in the spectrum to the right" exactly by the amount it itself is shifted from the origin. Vice-versa analogy holds for left shifted impulses. This is the mathematical background of this if you needed one. And yes, those results look correct to me. Cheers!


Wavelet transform in control systems In control systems, the Laplace transform is often used to analyze the stability and the performance of LTI system. For instance, the LTI system is stable if and only if the transfer function, which is the quotient between the Laplace transform of the output of the system and the Laplace transform of the input, has all of its poles in the left half complex plane.

Does wavelet transforms also found applications in analysis or design of control systems?

In the paper Haar-Based Stability Analysis of LPV Systems, the Haar wavelet transform theory have been used to design linear matrix inequalities (LMIs) to analyze the stability of LPV systems.

As the resolution level of the Haar wavelet increases, the number of variables and rows of the designed LMI increases and the feasibility of this LMI becomes a less conservative condition for the stability of the LPV system. Although the stability of LPV systems using LMIs that become less conservatives with the increase of variables and/or rows have already been proposed in the literature, the Haar-based approach can handle a larger class of parametric dependencies as well as non-convex parametric domains.


Barker Sequence for QPSK I'm trying to implement a QPSK burst modulation and demodulation through a channel and I'm now considering the synchronization part.

I've read that Barker sequence were excellent option for this as their autocorellation looks like a dirac.

So far I know that they're composed of 1 and -1 and that the longest found is 13 characters long.

My question is the following:

Is there some adaptation of a Barker sequence for QPSK to make use of
the complex representation or do you simply send a classic sequence
using only 1 and -1 while you could use also i and -i?

Polyphase Barker sequences have been studied (Google throws up multiple hits), but from a practical viewpoint, it is usually simpler to use binary Barker sequences on the I and Q channels.


A question about the Commutativity Property in discrete signals block diagrams According to the Commutativity Property of block diagrams, we can place b0 at the entrance or at the exit of the block diagram, it wouldn't matter.

My question is, if we placed it at the entrance, we should have the following difference equation:

$$y[n]=b_0x[n]-a_1y[n-1]$$

But if we placed it at the exit, we should have this difference equation:

$$y[n]=b_0x[n]-a_1b_0y[n-1]$$

Which don't correspond to the same transfer function, can anybody explain to me how it won't matter where to put b0 and that both schematics would lead to the same transfer function.

If you apply the multiplication with $b_0$ at the output, then that multiplier is after the feedback loop. Consequently, if the output (after the multiplier) is $y[n]$ then the input to the delay element is $y[n]/b_0$. The corresponding difference equation is

$$y[n]=b_0\left(x[n]-\frac{a_1}{b_0}y[n-1]\right)=b_0x[n]-a_1y[n-1]\tag{1}$$

which is identical to the one with the multiplier at the input, as expected.


Strange noise deviation formula I'm having troubles trying to use one relatively simple noise estimation method, here it is:A Fast Method For Image Noise Estimation Using
Laplacian Operator and Adaptive Edge Detection

My problem is with their proposed noise deviation formula (page $2$, ($3$)), I can't seem to comprehend it. What I want to know am I right that this $N$ is amount of pixels in image?

And also I'm having troubles implementing first part of the formula since its result is float undeflow. Maybe I can just use regular pixel_value - pixel_mean formula? Will it produce same results?

Altrough images with massive noise will have mean that actually pretty much depends on the noise amount and thus estimation will be considerably lower than the real value.

Where can I read about this formula and how to use, or how it's even named?

Maybe I've implemented it wrongly, here's my code that does it:
noise deviation implementation (C++ CImg)
It gives me very big numbers, like for noisy image of bears I got something like 1830152.

First of all, what I'm supposed to get from that formula? Isn't it supposed to be some kind of 0.xx coefficient or maybe percent of noise in image?

I can answer the first part of your question: Where formula 2(3) comes from.

Since the noise model is Laplacian, you are looking at an estimate of the standard deviation of a Laplace distribution. Compare the moments equation on this Wikipedia page: http://en.wikipedia.org/wiki/Laplace_distribution#Moments

As for what N is, it looks like they define it as a matrix operator in 2(2). So the sum in 2(3) is taken over pixels of an image generated by convolving the original image I(x,y) with the Laplacian operator N (which is denoted by $I(x,y)*N$, where $*$ is the convolution operator).


Measure similarity between quantized grayscale images I want to measure the similarity between quantized grayscale images (8 levels) that mostly depict organic patterns found in nature, say pigmentation patterns in animals (zebra pattern, leopard spots,..). I cannot use some kind of trained model, and I need a continuous similarity value. Feature scale tolerance is not important, but differences in rotation and density of the features could be large. Ideally two different flat textures of zebra patterns should have a similarity of 1.0 and say zebra-tiger a similarity of 0.7. I understand how impossible that sounds, but I hope someone could point me to relevant current research or methods somewhat suitable to what I'm looking for.

EDIT: More info about the desired usage.

Examples of images undergone grayscale conversion and quantization:
(512x512 pixels, 8 gray levels). I do not own these images.

These are then used as templates for a procedural pattern generation technique I am developing. Images are randomly generated, generation settings that produce high similarity are kept to be further tweaked. The procedurally generated images could be anything from random pixels to intricate shapes, slowly converging at the image of highest possible similarity.

This question is very broad, yet interesting. Your way to go are co-occurence matrices.

Co-occurence matrices can be calculated for each image and provide further indicators which specifically allow identification of different textures.

Note that there does not exist a single co-occurence matrix for an image, but several ones, since the co-occurence matrix is defined for specific directions of occurence. Further, you can vary the granularity of the co-occurence matrix. For each co-occurence matrix you can derive several features which enable you to distinguish the textures.

This example is very simple, yet should enable you to understand:

Image you have black-white picture of only a single line of pixels such as:

[0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0]

Now there are there are four possible co-occurences of pixel neighbours:

1) [0,0]
2) [0,1]
3) [1,0]
4) [1,1]

For each these we will now count the times that you can them within the strings:

1) [0,0] ---> 4
2) [0,1] ---> 3
3) [1,0] ---> 3
4) [1,1] ---> 0

This will give us our co-occurence matrix:

[[4,3],[3,0]] which is specific for each pattern.

Based on this co-occurence matrix, you are able to calculate features, such as Contrast, Correlation, Energy, Homogenity and many more.

Matlab provides you with tools to get the co-occurence matrix. Further, wikipedia provides an introduction.

Once you implemented this and encountered your first problems, I am looking forward to your next question. Good luck.


Extract frequencies from PCM signal I have an PCM (pulse-code modulation) signal. If I want to work out the frequency of a sample $ x $ long, I count the number of peaks/troughs of the sample, and divide by $ x $. The equivent dB of frequencies $ 20Hz > 20 000 Hz $ is then just the equivelent dB of the signal as a whole

However, I would like to create a spectrum analyzer of the type shown below. This requires extracting the equivelent decibels of different frequencies in the signal. So I need to split the signal into different signals of different frequencies, so I have signals of $ 20Hz > 5000Hz $, $ 5000Hz > 10000Hz $,  $ 10000Hz > 20000 Hz $, ect.

How can I extract signals for the different frequencies from one master signal?

It is not clear if you want to visualize the signal or process/transform it ("graphic equalizer" is an effect/transformation adjusting the signal level independently in various frequency bands, not a visualization - a more correct term for the representation you posted an image of is "spectrum analyzer").

It is not clear if you are working in the analog domain or the digital domain, and what kind of technology you want to employ (analog electronics? computer code running on a DSP chip or microcontroller? on a desktop computer?).

Are your signals audio signals? In which case you have to be aware that a meaningful representation should display the 20 Hz .. 20kHz range (range of human hearing), and that it should be subdivided logarithmically.

To get this kind of visualization, there are two solutions:

A bunch of band-pass filters (and a low-pass + a high-pass for the lower and higher channels). You then measure the (log-) energy in each band to get the "height" of each bar. This can be done in the analog domain (a bunch of op-amps for the filter, a rectifier to measure the signal level, and a log-amplifier to get the energy on a log scale; or with dedicated ICs like the MSGEQ7); or in the digital domain with digital filters.
A FFT + some grouping of adjacent bins (which is more or less equivalent to computing the filter bank in the frequency domain). Note that this can only be used for digital signals.


How to Prove a System Is Invertible? what i know is that for a system to be invertibel it should be one-one ,
but I am confused that if i am given a transfer function of a LTI system how can I prove or verify if it is invertible.

Another query is that can a system has more than one inverse systems , sometime before I used to think that inverse pairs are unique but just solving an example of finding inverse of differentiation as LTI system i came to result that it has got two inverses

1.u(t)

2.-u(-t)

In general LTI System is invertible if it has neither zeros nor poles in the Fourier Domain (Its spectrum).
The way to prove it is to calculate the Fourier Transform of its Impulse Response.

The intuition is simple, if it has no zeros in the frequency domain one could calculate its inverse (Element wise inverse) in the frequency domain.

Few remarks for the practice world:

If for a certain system we know the input signals are limited to certain band (Or zones) in frequency domain it is enough the LTI system has no zeros in this zone only.
In practice we have noise hence even values which are not mathematically zero but close to 0 (Or even low relative to the gain of the LTI system, Say ~-60 [dB]) are making the system not invertible.

Look for many question in this community about Deconvolution and Inverse Problems.


Is fft2 in MATLAB unitary? Some differences happen I meet a problem when implementing fft2 in MATLAB.

The question is I try to simulate the realistic measurements $Y = |FCXF^H|^2$ - the intensity of Fourier domain of object $X$, where $F$ denotes the Fourier transform matrix $C$ denotes $\{+1, -1\}$ matrix.

But I want to replace the $C$ to realistic $D = \{0,1\}$ with $2D - 1 = C$. Therefore, after some simple calculations, the measuremnts $Y$ becomes $Y_{real} = 2(|FDXF^H|^2 + |F\bar{D}XF^H|^2) - |FXF^H|^2$, where $\bar{D}$ denotes the inversion of $D$. The calculations are from a paper.

For original $Y$ and $C$ the algorithm could work. Therefore, I try to verify with $Y_{real}$. Strange things happen, although the intensity of Fourier domain of $Y_{real}$ looks similar with $Y$. They are totally different!!!

I am very confused, the equations are right, but why they are different? I think the problem is from fft2 function. Could anyone tell me the principle of fft2? Thanks in advance!!!

Here is the code:

This implements the original algorithm with $C$:

Y = abs(fftshift(fft2((2*D-ones(n1,n2)).*x))).^2;

This implements the equation - realisitic measurement:

Y_real = 2*(  abs(fftshift(fft2(D.*x))).^2 + abs(fftshift(fft2( ( ones(n1,n2) - D ).*x ))).^2 ) - abs(fftshift(fft2(x))).^2;

The standard (conventional) definition of DFT (1D or 2D) is not unitary.

See for example the 1D standard (conventional) DFT pair as:

$$ X[k] = \sum_{n=0}^{N-1} x[n] e^{-j \frac{2\pi}{N} n k } $$

and

$$ x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j \frac{2\pi}{N} k n} $$

DFT is not unitary due to the fact that the forward and backward transforms are not symmetric (the scale $\frac{1}{N}$)

Hence MATLAB's fft() and fft2() functions will not provide unitary transforms; i.e., they do not preserve energy.

However, you can define unitary version of DFT by distributing the scale as:

$$ X[k] = \frac{1}{\sqrt{N}} \sum_{n=0}^{N-1} x[n] e^{-j \frac{2\pi}{N} n k } $$

and

$$ x[n] = \frac{1}{\sqrt{N}} \sum_{k=0}^{N-1} X[k] e^{j \frac{2\pi}{N} k n} $$

Now this pair is unitary.


Good metric for qualitatively comparing image patches I am trying to "match" little square patches in an image. At first glance, it seems reasonable to simply do a Euclidean distance style comparison of two of these arrays to get a "similarity" measure. This works fine in many cases (the "best" patch (lowest value) according to this metric looks very much like the query patch). However, there are many cases in which this produces a very bad match. For example, take these two patch pairs:

Two patches of a brick wall, score 134 (this is the sum of the components of the average absolute pixel difference):

One patch of a brick wall, one patch of grass, score 123!

To a human, "clearly" the grass does not match the brick, but this metric says otherwise. The problem is just in the local statistical variation.

If I use something like a histogram comparison, I completely lose all spatial information - e.g. if a patch is grass on the top and brick on the bottom, it would match exactly to a patch with grass on the bottom and brick on the top (again, another "obviously wrong" match).

Is there a metric which somehow combines both of these ideas into a reasonable value that will evaluate to "similar" for Pair 1 above, but not also be similar for my example of a patch and its vertical mirror?

Any suggestions are appreciated!

To start with, 2 images and 3 measurements are not exactly enough to ascertain any kind of statistical model which in terms can be used to define the optimal comparison metric.

I think you could start to have a look at texture recognition papers for methods and clues. It is an active field.

For what it is worth, I ran a couple of perceptual hashing functions (DCT and Random Projection based) tests and did a little trial with the so called SIFT descriptor. These functions can separate inter and intra class distances, although from 3 images it is impossible to conclude anything.

Code on github.


Moving Average for Notch filtering I have a periodic signal(ECG) with period of ~1 seconds. It does not have features that are shorter than 0.04 seconds. For removal of 60Hz, I thought instead of implementing a notch filter, doing a simple moving average of 1/60~0.016 seconds.(Averaging over a larger window might destroy signal information, and also it's not needed) This is much simpler than a notch filter which can be of several orders. Am I missing something?

No you are not missing anything. A moving average with a period of $T = 1/60$ seconds will indeed have a notch at 60 Hz, since the frequency response is a Sinc function with the first null at 1/T.  (So it is both a low pass filter with an magnitude envelope that goes down as 1/f combined with a notches at integers of 1/T not including 0).

This is because a moving average has an impulse response that is a rectangular function of duration T, and the Fourier Transform of such a function is a Sinc (and the Fourier Transform of the Impulse Response is the Frequency Response). Digital Filters are similarly formed when the impulse response is a sampled rectangular function and thus their Fourier Transform will approach a Sinc function (what I like to refer as an "aliased" Sinc function similar to other aliasing in discrete time systems). If we position the sampling rate at integer multiples of 60 Hz, starting at 120 Hz and higher, the null will be at 60 Hz and it's higher harmonics if the digital impulse response (formally called the unit sample response) is uniform samples for duration T, where T is 1/60 Hz (which is a moving average specifically when the weighting of each sample is 1/N for N samples, otherwise it is a scaled moving average). Therefore the minimum sampling rate required to do this digitally would be 120 Hz, otherwise 180 Hz, 240 Hz, and higher multiples of 60 Hz can also be used).  Note that I said this would approach a Sinc function; notice the plots below of the actual frequency responses for each choice of sampling rate when we position the null at 60 Hz and do the moving average, given by the actual digital transfer function according to:

$$H(z) = \frac{1}{N}\sum_{n=0}^{N-1}z^{-n}$$

With the frequency responses below determined by using the unit circle for z ($z= e^{j\omega}) $ where $\omega$ is the normalized angular frequency from 0 to $2\pi$ corresponding to 0 to the sampling rate. The actual Matlab/ octave command that can be used to readily determine this is as follows (and Python SciPy has a similar command):

freqz([ones(N,1), 1/N]);

Where N for each of the cases below is 2, 3, 4, and 10.

Also what you could alternatively consider which is even simpler is a "comb filter" structure which is simply the addition of the input with a delayed version of itself where the delay is equal to T/2.

This will have the frequency response according to:

$$H(\omega) = 1 + e^{-j\omega T/2}$$

Since the Fourier Transform of a delay is given as

$$x(t) = \delta(t-T) <==> x(\omega) = e^{-j\omega T}$$

With the magnitude plotted below:

This is equivalent to the digital filter with a delay of one sample when the sampling rate is 120 Hz, and a transfer function given as:

$$H(z) = 1+z^{-1}$$

with as stated previously the frequency response when z is the unit circle as

$$H(e^{j\omega}) = 1 + e^{-j\omega}$$

with $\omega$ going from 0 to $2\pi$ in this case representing the normalized angular freq from 0 to the sampling rate.

Like the moving average case above, this could also be implemented at any sampling frequency that is an integer multiple of 60 Hz, starting at the minimum 120 Hz rate.

If you are concerned about the "frequency droop" in your passband area for either of the above two solutions, you could consider using this very simple three tap FIR compensator that I described at this link to help flatten the passband:

how to make CIC compensation filter

In all the cases above, the sampling rate must be an integer multiple of the notch frequency desired. If that is too restrictive, here is a digital notch filter that I would typically implement when a very tight notch is required with a flat passband, to compare to the approaches above (and does not have the same limitations on sampling rate):

Transfer function of second order notch filter

Side note: The relatively high sidelobes, or stated another way, the slow frequency roll-off of 1/f in magnitude, is why a moving average is often a poor choice as an FIR filter for low pass filtering (or estimating the DC value which IS the average) UNLESS you are filtering what is known to be white noise. For the case of white noise it is one of the best choices. Why? Because the moving average does have the narrowest main lobe for the given number of taps compared to any other FIR filter, and if the noise is known to be white, then the 1/f magnitude roll-off will be sufficient to suppress the higher frequency components.


Real and complex low pass equivalent of band-pass signal A real band-pass signal $s(t)$ with center frequency $f_c$ can always be represented as

$$ s(t) = A(t) \cos(2\pi f_c t) $$

where $A(t)$ is real and can be viewed as amplitude.

Also

$$ s(t) = A(t) \cos(2\pi f_ct) = \Re\left\{(A(t) e^{j2\pi f_ct})\right\}$$

Then we consider

$$s_l(t)= A(t) e^{j2\pi f_ct}$$

to be the corresponding complex band-pass signal of $s(t)$.

My question is how can I find the real low pass equivalent of $s(t)$ and the complex low pass equivalent of $s_l(t)$?

It's not true that any real-valued band pass signal can be written as shown in your first equation:

$$s(t)=A(t)\cos(\omega_c t)\tag{1}$$

The most general representation of a real-valued band pass signal is given by

$$s(t)=\Re\{a(t)e^{j\omega_c t}\}\tag{2}$$

where $a(t)$ is a complex-valued function, called the complex envelop of $s(t)$. It is the complex baseband (low pass) representation of $s(t)$. With $a(t)=a_I(t)+ja_Q(t)$, (2) can be written as

$$s(t)=a_I(t)\cos(\omega_c t)-a_Q(t)\sin(\omega_c t)\tag{3}$$

where $a_I(t)$ and $a_Q(t)$ are called the in-phase and the quadrature components, respectively. Both $a_I(t)$ and $a_Q(t)$ are real-valued baseband (low pass) signals.

A third, equivalent, representation of $s(t)$ is obtained by writing the complex envelop as

$$a(t)=b(t)e^{j\phi(t)}\tag{4}$$

where $b(t)$ is the (real-valued) envelop, and $\phi(t)$ is the phase. Using (4) and (2), the band pass signal $s(t)$ can be written as

$$s(t)=b(t)\cos(\omega_ct+\phi(t))\tag{5}$$

Equations (2), (3) and (5) are three equivalent representations of a general real-valued band pass signal. If you compare representation (5) to Equation (1), you see that Equation (1) is a special case for $\phi(t)=2\pi k$, $k=0,1,\ldots$, i.e. for a constant phase which is zero (or a multiple of $2\pi$). From Equation (3), this is equivalent with the condition that the quadrature component $a_Q(t)$ is zero.

Note that the complex-valued band pass signal

$$\tilde{s}(t)=a(t)e^{j\omega_c t}\tag{6}$$

is the analytic signal mentioned in a comment. It has no negative frequency components and it can be obtained from $s(t)$ by computing

$$\tilde{s}(t)=s(t)+j\mathcal{H}\{s(t)\}\tag{7}$$

where $\mathcal{H}\{\cdot\}$ denotes the Hilbert transform.


Counting the number of groups of 1s in a boolean map of numpy.array I am right now dealing with some image processing in Python via PIL (Python Image Library). My main aim is counting the number of colored cells in an immunohistochemistry image. I know that there are relevant programs, libraries, functions and tutorials about it, and I checked almost all of them. My major aim is writing the code manually from scratch, as much as possible. Hence I am trying to avoid using lots of exterior libraries and functions. I have written the most of the program. So here is what's going on step by step:

Program takes in the image file:

And processes it for the red cells (basically, it turns off the RGB values below a certain threshold for red):

And creates the boolean map of it, (gonna paste a part of it since it is big) which  basically just puts 1 wherever it encounters with a red pixel in the processed second image above.

22222222222222222222222222222222222222222
20000000111111110000000000000000000000002
20000000111111110000000000000000000000002
20000000111111110000000000000000000000002
20000000011111100000000000000000001100002
20000000001111100000000000000000011111002
20000000000110000000000000000000011111002
20000000000000000000000000000000111111002
20000000000000000000000000000000111111102
20000000000000000000000000000001111111102
20000000000000000000000000000001111111102
20000000000000000000000000000000111111002
20000000000000000000000000000000010000002
20000000000000000000000000000000000000002
22222222222222222222222222222222222222222

I intentionally generated that frame like thing on the borders with 2s to help me with counting the number of groups of 1s in that boolean map.

My question to you guys is, how come I can efficiently count the number of cells (groups of 1s) in that kind of boolean map? I have found http://en.wikipedia.org/wiki/Connected-component_labeling which look extremely related and similar, but as far as I see, it is at the pixel level. Mine is at the boolean level. Just 1s and 0s.

Thanks a lot.

You can use ndimage.label, which is a nice way to do this. It returns a new array, with every feature having a unique value, and the number of features. You can also specify a connection element.

import scipy
from scipy import ndimage
import matplotlib.pyplot as plt

#flatten to make greyscale, using your second red-black image as input.
im = scipy.misc.imread('blobs.jpg',flatten=1)
#smooth and threshold as image has compression artifacts (jpg)
im = ndimage.gaussian_filter(im, 2)
im[im<10]=0
blobs, number_of_blobs = ndimage.label(im)
print 'Number of blobls:', number_of_blobs

plt.imshow(blobs)
plt.show()

#Output is:
Number of blobls: 30