Phase reference of a periodic signal Assume an arbitrary (discrete) signal that is periodic and known over a whole period. I need a way to select a characteristic point along the signal such that I can always retrieve it even when the signal is phase shifted.

Candidates could be zero crossings or extrema or inflections... provided they exist and are unique or can be discriminated from others.

For example, the position of the maximum would work well for a triangle wave, a little less accurately for sinusoid and would raise an ambiguity for $\sin x+\frac13\sin 3x$ as this one has two maxima of equal height.

The method should be robust to noise and work for relatively smooth signals.

Do you know any general feature/formula having the desired properties ?

Compute a DFT of your signal over a period of $T$ samples, extract the phase $\phi$ (in the $[0, 2\pi[$ range) of the first harmonic, and use the point $\frac{2\pi-\phi}{2 \pi}T$ as your landmark.

This will not correspond to a particular characteristic landmark in the time-domain waveform; but the procedure to retrieve it is by design invariant to phase shifts, distortion adding higher harmonics, and to a large extent noise.

If the waveform of your signal is known in advance you can also use it as a template, and select as a landmark the lag, in samples, at which the cross-correlation between your observed period and your template is maximal.

Python code here:

import numpy
import pylab

def shift(s, k):
  n = s.shape[0]
  return s[(numpy.arange(n) + k) % n]

T = 256
t = numpy.arange(0.0, T) / T

original = numpy.sin(t * 2 * numpy.pi) + numpy.sin(3 * t * 2 * numpy.pi) / 3
time_shifted = shift(original, 40)
distorted = shift(original + original ** 3, 200)
noisy = shift(original + numpy.random.randn(1, T).ravel() * 0.2, 180)

signals = [original, time_shifted, distorted, noisy]

for i, s in enumerate(signals):
  pylab.subplot(411 + i)
  pylab.plot(s)

  # Phase of first harmonic method.
  first_harmonic_coefficient = numpy.fft.rfft(s)[1]
  angle_2pi = -numpy.angle(first_harmonic_coefficient) + 2 * numpy.pi
  angle_2pi = angle_2pi % (2 * numpy.pi)
  landmark = numpy.round(angle_2pi / (2 * numpy.pi) * T)
  pylab.plot(landmark, 0, 'r+')

  # Cross-correlation method.
  xcorr = numpy.fft.irfft(numpy.fft.rfft(s) * numpy.fft.rfft(original))
  landmark = numpy.argmax(xcorr)
  pylab.plot(landmark, 0, 'g*')

pylab.savefig('plot.png')


Calculate amplitude of noise added Sine Gaussian signal from frequency domain How to calculate the time domain amplitude of noise added Sine Gaussian signal from corresponding frequency domain. Here time domain amplitude is 2. How to get this from frequency domain.

N = 4096;
Samplingfrequency= 1024;
timelength=N/Samplingfrequency;

% Generate random
Noise=rand(1,N)*.1;

t    = linspace(-timelength/2,timelength/2,N);
Freq=(1/timelength)*[0:N/2-1 -N/2:-1];

signalFrequency = 64 ;
amplitude=2;
Q=10;
tau  = Q /(2 * pi * signalFrequency);

%Sine Gaussian
Gaussian=((amplitude*exp(-1 .* ((t)/tau).^2) .* cos(2 * pi * signalFrequency .* t))+Noise);

%Time Domain
subplot(2,1 ,1)
plot(t,Gaussian);

%Frequency Domain
FFTGaussian=abs(fft(Gaussian));

subplot(2,1 ,2)
% One sided frequency domain plot
plot(Freq(2:N/2),FFTGaussian(2:N/2));

Random noise doesn't have an amplitude, it has statistical properties such as standard deviation and variance.  Variance (standard dev. squared) is interpreted as the AC power of the noise.  In the freq domain you can estimate the total power of a signal and the total 'power' (variance) of the noise.  The ratio of the signal power over the noise power is your signal-to-noise ratio.  Multiplying the log-base-ten of that ratio by ten gives you the signal-to-noise ratio measured in dB.  [-Rick-]


A signal containing no pure sinusoidal components I'm trying to decipher what exactly this statement means.  So if I have a signal which has no pure sinusoidal components, does that mean I can't decompose my signal completely into sine-waves which may be amplitude/phase/period shifted?  And so I would conclude that there is also like a saw-tooth component or a square-function component mixed up in the signal as well?

Saying that the signal contains NO pure sinusoidal component makes me feel like its decomposition should not include ANY amplitude/phase/period shifted sine-waves.. but that seems a little extreme.

Could someone clarify this for me?  Thanks.

It contains no spikes in its frequency spectrum. (What is the Fourier transform of a sinusoid?)


Intensity change after filtering this image? I need to filter the image raster below with the following mask:

 0 1
-1 0

and the Image Raster:

0 0 0 0 0 1
0 0 0 0 1 1
0 0 0 1 1 1
0 0 1 1 1 1
1 1 1 1 1 1

My results are:

 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
-1 0 0 0 0

In what direction does this mask detect intensity change? What does the sign of the output signal mean?

For binary images, the filter kernel you have is looking for a $0 \rightarrow 1$ edge in the North-East direction.  Your image raster only has one place where such an edge exists (in the bottom left / South-West corner).

The sign of the output means that the edge is a $1 \rightarrow 0$ edge transition rather than a $0 \rightarrow 1$ transition.


Scaling down a raster image I have a raster image (A) of size 20x20 and I want to scale it down to a size of 10x10 (B).

Naturally in the resulting picture B one pixel will represent 4 pixels shaped 2x2 from A.

Is it possible to give a canonical answer to how the pixel values of B have to be calculated if no neighbouring pixels (surrounding the 2x2 subset) in A should be taken into account? Mean, maximum or something else?

Bilinear is the most widely used method. The nearest neighbor down-sampling algorithm is the fastest but least accurate.

Note that when you are trying to down-scale an image by half, the bilinear sampling algorithm becomes the average (mean as you mentioned):

x1[i/2][j/2] = (x[i][j] + x[i+1][j] + x[i][j+1] + x[i+1][j+1]) / 4;


Difference Equations, why M<=N for causality? in my notes for DSP they have the difference equation general form as:

y(k) + a1*y(k-1) + a2*y(k-2) + ... + an*y(k-N) = b0*x(k) + b1*x(k-1) + ... + bm*x(k-M)

with the claim that for the output to not depend on future values of the input then M<=N.

Why is this the case?

In fact, if I'm not missing something here I think I can provide a counterexample:

y(k) = x(k-1)

which has N = 0 and M = 1.
But this is just a unit delay, the output is just the previous value of the input.

What's going on here?

You're right, it's simply not true that $M\le N$ is necessary for causality. The difference equation in your question can always be implemented by a causal system.

However, note that the difference equation is not uniquely related to a causal system. A simple example is

$$y[n]=ay[n-1]+bx[n]\tag{1}$$

This can obviously be implemented by a causal system. But for $a\neq 0$ you can also rewrite (1) as

$$y[n-1]=\frac{1}{a}(y[n]-bx[n])\tag{2}$$

which suggests an anti-causal system, even though (1) and (2) are completely equivalent. This difference is reflected by the transfer functions ($\mathcal{Z}$-transforms) of the corresponding systems. For a causal system the region of convergence (ROC) is outside a circle enclosing all poles, whereas for an anti-causal system the ROC is inside a circle outside of which all poles are located.


Far field spherical waves plane waves All em radiation propagates spherically and at a far off place a small enough distance, when observed appears like a plane like our experience on the earth.  So I tend to think that whether a wave on an aperture can be considered plane or not depends on how far the aperture is from the origin and how much is the aperture.  So it should be only dependent on the angle made at the centre by the aperture i.e. apertures making more than a certain angle at the centre can be considered spherical, and otherwise can be considered plane.  This is my wrong understanding.    But in radar theory, this is not the criterion for considering plane and far field.  A 2D^2/lambda formula is used.  How is this formula relevant?  How lambda is entering into the formula?

Rajeev Bansal in the appendix to this article explains the $2D/\lambda$ rule:

How does one arrive at the various criteria for the far-field zone?
Basically, the criteria are guidelines to the boundaries where the
fields start to approximate the “ideal” assumed characteristics. First
of all, only the radiation ($1/r$) terms remain significant; higher
order terms fade away. Second, in the far-field zone, the angular
field distribution becomes independent of the distance [5]. Third,
only transverse field components remain, and the ratio of the electric
and magnetic field components approaches the free-space impedance, 377
ohms [6 ,7]. Finally, for a receiving antenna, the incoming wave-front
is nearly planar across the aperture. In fact, the $r = 2 D^2/\lambda$
formula corresponds to a phase error (due to the curvature of the
actual spherical wave-front) of no more than $\pm 22.5$ degrees across the
aperture, as compared with the ideal plane wave-front [6, 8].
References:

IEEE Standard Dictionary of Electrical and Electronics Terms, IEEE standard 100-1984, 1984.

C. Paul, K. Whites, and S. Nasar, Introduction to Electromagnetic Fields, 3rd ed., McGraw-Hill, 1998.

Reference Data for Radio Engineers, 6th ed., SAMS, 1981.

G. Smith, An Introduction to Classical Electromagnetic Radiation, Cambridge, 1997.


Comparing two Accelerometer signals in time domain As a research to answer my question, I've been reading several posts on this stack with a common thread of "comparing signals". The methods recurrently suggested were to use DTW, correlation or MSE etc.
The two signals I'm trying to compare are generated by piezoelectric single-axis accelerometer and they are periodic.
Out of 4 variables playing a role in signal amplitudes, two of the variables are different(A,B resp) and the rest are same(C,D resp).
My aim is to measure quantitatively how different the two signals are on the basis of variable B. At the moment, a qualitative way to go about this, according to my supervisors, is to look which among the two signals have the higher peaks.
I want to get a quantitative measure to compare the two signal's peaks or maybe the whole dataset which consists of 20,000 datapoints for a period of 10 sec. Are peaks/amplitudes the right way to go about this?

I suggest using a correlation metric as a comparison. Remove any mean in the signals and after subtracting the mean, compare the signals using the Pearson Correlation Coefficient to have a single data point that scales between +1 and -1 with regards to their direct correlation assuming any time offset factors have been removed, and even better as a first evaluation use the cross-correlation function which would provide the above (without normalization) for every offset in time, in case their are delays between the results.


Matched filter for "amplitude SNR" vs power SNR I am working on an application for which a matched filter seems like the right concept.  In the derivation for the matched filter which I went through (here), they define the SNR as the ratio of signal power to noise power as the basis for the derivation of the expression for the matched filter.  I think that for my application, the "amplitude" SNR would be a better target for optimization.  However, it looks like I reach a mathematical dead-end that is avoided by the use of the power SNR.  Below is my work which roughly follows the derivation from wiki:

Is there a manipulation I am missing that could get me past this obstacle or does someone know if this result has been solved for already?  I have not found it anywhere.

Power is proportional to the square of amplitude, so there's really no difference (since SNR is always a positive number) when maximizing: A filter that maximizes "amplitude SNR" maximizes "power SNR".

So, if your application requires you to maximize the amplitude ratio between signal and noise (be careful, amplitudes can be negative or complex, depending on how you define them!), then that's the same as requiring you to maximize the power ratio.

I.e.

$$h_\text{opt} = \arg\max_h E\left\{\left(h^H (vv^H)h\right)^{\frac12}\right\} = \arg\max_h E\left\{h^H (vv^H)h\right\} $$

since the quadratic form $h^H (vv^H)h$ (which is the power of $y$) can never be negative, and that makes your square root both unambiguously invertible and monotonous.


Corner detection using Chris Harris & Mike Stephens I am not able to understand the formula,

What is $W$ (window) and intensity in the formula mean,

I found this formula in opencv doc http://docs.opencv.org/trunk/doc/py_tutorials/py_feature2d/py_features_harris/py_features_harris.html

Harris Corner detector tries to quantify the local intensity changes at all the directions for each pixel. The figure below illustrates the basic idea clearly:

So $I(x+u,y+v)$ indicates the pixel intensities of all the neighborhood pixels around $(x,y)$.
The window function is applied for feature localization. For most often used Gaussian function, the choice of sigma value corresponds to different feature width in the image.


Convolution of two exponential signals I'm a little new to signal processing and I'm trying to wrap my head around convolutions.

I know the definition of convolution for a continuous signal is

$$y(t) = x(t) * h(t) = \int_{-\infty}^{\infty}{x(\tau)h(t-\tau) \, \mathrm{d}\tau}$$

Let's say for example that

$$h(t)=6e^{-t}u(t)$$

and

$$x(t)=e^{-4t}u(t)$$

you end up getting that

$$y(t)=\int_{-\infty}^{\infty}{6e^{-3\tau-t}}d\tau$$

which is divergent. I watched a few YouTube videos about it and they always explain how to do convolutions graphically with two rectangles, does anybody mind explaining how to do it with two exponential functions, or point me to somewhere that does?

Just in case if you look for a complete solution based on integration:

$$y(t)=\int_{-\infty}^{+\infty}e^{-4\tau}u(\tau)\times 6e^{-(t-\tau)}u(t-\tau)d\tau$$
Consider these facts:

Both step functions need to be equal to $1$ to have both terms (otherwise the corresponding term, and consequently $y(t)$, become zero). This implies:

$\tau>0$
$t-\tau>0\Rightarrow \tau<t$
From the two above, $t>0$. Otherwise the integral is zero.

We need $1$ and $2$ for the boundaries of integral. They say $\tau \in (0,t)$. Also we need $3$ to know where the output exists. It simply says the integral should be multiplied by $u(t)$. Now that we applied all the given information from the two step functions, we remove them and the integral becomes:
$$\begin{align}
y(t)&=\left(\int_{0}^{t}e^{-4\tau}\times 6e^{-(t-\tau)}d\tau\right)u(t)\\[10pt]
&=\left(6e^{-t}\int_{0}^{t}e^{-3\tau}d\tau\right)u(t)\\[10pt]
&=\left(6e^{-t}\left(\frac{e^{-3\tau}}{-3}\biggr\vert^t_0\right)\right)u(t)\\[10pt]
&=\left(-2e^{-t}\left(e^{-3t}-1\right)\right)u(t)\\[10pt]
&=2(e^{-t}-e^{-4t})u(t)
\end{align}$$


Real time fft/ifft processing when manipulating a live stream of audio data in the frequency domain, I've heard you need to do fft / ifft in small chunks, say 256 samples at a time.

I've also heard you need to somehow employ a window function but I'm not sure the details on that. Is that to prevent artifacts from popping up due to the circular nature of fft?

Is there anything else I should know that I haven't mentioned above if planning to do this? I also eventually want to figure out how to make windows work with convolution reverb, but assume the answer to that is big enough for its own question.

Thanks for any info!!

Analysis is separate from processing.  Windowing helps remove certain artifacts during spectrum analysis, but since it is numerically lossy, isn't as useful when actually processing the audio for output.  Instead, sufficient zero-padding (without any non-rectangular windowing) and the use of overlap-add or overlap-save is used for any fast convolution or frequency domain filtering.  This prevents (actually just zeros) circular convolution artifacts.

Convolution reverb is already answered elsewhere, but in general can be done by breaking up the reverb impulse response into a whole series of overlap-save processing steps and summations (e.g. a 256k long reverb might be done using a series of 1,000 256+padding length FFT convolutions, or a mix of shorter and longer ones, with a future use buffer to hold most of those results till needed.)


What's the difference between nonstationary and time varying? I am confused by the following statement from the paper Beamforming: A Versatile Approach to Spatial Filtering by Barry D. and Kevin M.:
There are two basic adaptive approach: 1) block adaption ..., 2) continuous adaption.... If a nonstationary environment is anticipated, block adaption can be used, provided that the weights are recomputed periodically....Continuous adaption is usually preferred when statistics are time-varying or (for computational reasons) when the number of adaptive weights M is moderate to large (values of M>50 is not uncommon).
This statement gives the conditions when the two kinds of adaption methods should be used. But the condition NONSTATIONARY and TIME-VARYING, in my opinion, are the same. Can you point out in what situations they are different? Illustrative examples would be greatly appreciated.
Thank you.

My understanding :
Non-Stationarity is mentioned in statistical sense where channel coefficients are non-stationary random processes. You can look for Wide-Sense Stationary Processes and Strict-Sense Stationary Processes and their definition and properties.
Time-Varying is with regards to time-domain channel impulse response characteristics. Time-variance indicates that depending upon when any specific input is provided to the channel, the output will change because the channel is varying with time.
A particular channel model can be time-varying but wide-sense stationary. This would mean, depending upon when any specific input is provided to the channel, the output changes but the statistical properties of the output like correlation function, mean etc will follow wide-sense stationarity conditions.


Refinement equation and orthogonal wavelets basis I have a wavelet function space defined by,

$ \phi(x) = \sqrt(2) \sum_n h_{\phi}(n) \phi(2x-n) $ .

given the Fourier transform of the function is:

$\hat{\phi}(\omega) = \frac{1}{\sqrt(2)} \hat{h}_{\phi}(\omega) \cdot \hat{\phi}(\frac{\omega}{2})$

why does the sum of the dc component in jumps of $2\pi$ equals to 1, i.e.:

$\sum_{l=-\infty}^{\infty} |\hat{\phi}(\omega +2\pi l)|^2 = 1$

or even its necessarily constant. Either I am missing something trivial, or whatever, but I can't seem to agree with it.

thanks.

This follows from the orthonormality of shifted versions of $\phi(x)$:

$$\langle \phi(x),\phi(x+n)\rangle =\delta[n]\tag{1}$$

where $\delta[n]$ is the unit impulse. Taking the discrete-time Fourier transform (DTFT) of $(1)$ gives the equation in your question.

Note that the left-hand side of $(1)$ is the auto-correlation of $\phi(x)$ evaluated at integers $n$. The Fourier transform of the auto-correlation is $\left|\hat{\phi}(\omega)\right|^2$, where $\hat{\phi}(\omega)$ is the Fourier transform of $\phi(x)$. Sampling the auto-correlation results in the infinite sum of shifted spectra.

EDIT: The equation can be derived as follows: define the auto-correlation of $\phi(t)$ as

$$p(\tau)=\int_{-\infty}^{\infty}\phi(t)\phi(t+\tau)dt\tag{2}$$

Realizing that $p(\tau)=(\phi_{-} \star \phi)(\tau)$ (where $\phi_{-}(t)=\phi(-t)$), the Fourier transform of $p(\tau)$ is

$$P(\omega)=\left|\hat{\phi}(\omega)\right|^2\tag{3}$$

Sampling $p(\tau)$ at integer $n$ results in a periodic continuation of the spectrum:

$$\textrm{DTFT}\{p(n)\}=\sum_{k=-\infty}^{\infty}P(\omega+2\pi k)\tag{4}$$

With $p(n)=\langle \phi(x),\phi(x+n)\rangle$ and with $(3)$ we finally get

$$\textrm{DTFT}\left\{\langle \phi(x),\phi(x+n)\rangle\right\}=\sum_{k=-\infty}^{\infty}\left|\hat{\phi}(\omega+2\pi k)\right|^2\tag{5}$$


Finding dependence between two sparse signals I am new to signal processing and would really appreciate your help. happy to provide more context/details.

In my research, I have a system that has many sensors $X_1, \ldots, X_N$ and using this information, the system makes a prediction on a value, $A$. For simplicity, let's assume that $A$ is binary, and please keep in mind that this is a real system I am trying to improve so that means that I can tell whether $A$ was right or wrong (the system can be checked against "ground truth").

Now each sensor makes measurement sporadically (not necessarily randomly, but definitely not some obvious pattern) for reasons that are too complicated to go into right now. This means that $X_1$ for example can resemble this: $[\textrm{2, 5, 1, na, na, 1, na, na, 89, $\dots$}]$ where each index of this vector is a time step, and $\textrm{'na'}$ represents times when the sensor did not take a measurement. This means that the vector actually contains two sets of information: the value of each individual measurement, and when the measurements were done. There is a third value (but maybe that's not important for now) which is about how much effor each sensor put into collecting this information.

Now my problem is that I want to measure the dependence between $X$'s. i.e. given a specific vector of measurement from, say, $X_4$, what is the probability that $X_{10}$ was influenced by $X_4$. Specifically, I am looking for a number (could be a probability, or a correlation coefficient) that would say that $X_{10}$ was influenced by $X_4$, and I want this number for all pairs of $X$'s.

What I am trying to do is to factor in these influences into the system that is trying to make a decision. Since I know the ground truth, I can tell how good my "dependence measure" is.

I have zero background in signal analysis and I thought that you all might know something here. I think this is an interesting general problem and I would really appreciate some help.

If you have enough sample measurement data, you could reconstruct a time series for each sensor pair that eliminates any time step where either one or both of the sensor pairs are na.  Then you are left with a time series for that pair without any na's and can proceed with standard time series approaches for calculating correlation.  Then repeat those steps for each sensor pair.

For example:

X_1 = [2, 5, 1, na, na, 1, na, na, 89]
X_2 = [na, 4, 9, 8, na, na, na, 2, 9]

after reconstruction becomes:

X_1 = [5, 1, 89]
X_2 = [4, 9, 9]


Unable to understand SNR equation in Matlab Can someone please help me understand this equation of SNR given in MATLAB:

SNR_dB = 1:0.1:7
SNR = 10.^(SNR_dB/10)

To express power $P$ in $\mathrm{dB}$ one performs this operation:

$$P_{\mathrm{dB}}=10\log_{10}(P)$$

The inverse operation would then be

$$\frac{P_{\mathrm{dB}}}{10}=\log_{10}(P) \implies 10^{\frac{P_{\mathrm{dB}}}{10}}=P$$

And that's exactly what you are doing in that MATLAB equation (i.e. expressing the SNR in linear scale from the data in decibels).

Note that the .^ means that the operation is being carried out element-wise.


Can Sparse Fourier transform be used for sparse signal in other domain I've read that sparse fast Fourier transform can be used to compute the Fourier transform of a signal that is sparse in frequency domain much faster compared to FFT. My question is that can SFFT be used for signal that is not sparse in frequency domain but in some other domains?

No.

What did you expect? For sparse algorithm to work, the signal has to be sparse. If it's sparse after a integral transform, that doesn't mean it's sparse under a different transform.

Also note that according to this answer, the benefit of the sparse fft only apply to very sparse signals - a factor of at least about 2000 between the zero and non-zero bins.


General form of dilations and translations of wavelet function In some papers, the dilations and translations of a wavelet function is written as follows:

$$\psi_{j,k}=\frac {1}{\sqrt {2^j}} \psi\left({2^{-j}t-k}\right),\quad\text{where $j$ and $k$ are integers}$$

However, in some literature, the above formula is written in a different way:

$$\psi_{j,k}=\ {\sqrt {2^j}} \psi\left({2^{j}x-k}\right), \quad\text{where $j$ and $k$ are integers}$$

Are these two formula equivalent?

Yes, they are equivalent, since $j\in\mathbb{Z}$. Note that your first definition becomes equal to the second when $j$ is substituded by $-j$. It is just a convention, if you consider a positive $j$ to be a up- or downscaling.

So, when writing a paper or reading a book, it is always important to understand the definition of the transform used.


Why does knowing the impulse response allow you to determine the output for any LTI system? Going over some notes on LTI systems and I’m quite confused about how knowing the impulse response of a system allows us to recover the systems’ response to any signal. In the textbook screenshot below, it looks like they just replaced the Dirac deltas with $f$ when defining $y(n)$.

(Adapted from this answer on dsp.SE)

The reason that the impulse response (also called the unit pulse response for discrete-time systems) determines the output for arbitrary input $x$ to an LTI system is that

  The output of a linear time-invariant system in response to input $x$ is
  the sum of scaled and time-delayed versions of the
  impulse response.

This is because of the linear property of the system: if the response to input signal $x_i$ is $y_i$, then the response to input $x_1+x_2+ x_3 + \cdots$ is $y_1+y_2 + y_3 + \cdots$.  So, let's decompose the input signal $x$ into a sum of scaled
unit pulse signals. The system response to the unit pulse signal
$\cdots, ~0, ~0, ~1, ~0, ~0, \cdots$ is
the impulse response or pulse response
$$h[0], ~h[1], \cdots, ~h[n], \cdots$$
and so by the scaling property the single input value $x[0]$,
or, if you prefer
$$x[0](\cdots, ~0, ~0, ~1, ~0,~ 0, \cdots)
=  \cdots ~0, ~0, ~x[0], ~0, ~0, \cdots$$
creates a response
$$x[0]h[0], ~~x[0]h[1], \cdots, ~~x[0]h[n], \cdots$$

Similarly, the single input value $x[1]$ or creates
$$x[1](\cdots, ~0, ~0, ~0, ~1,~ 0, \cdots)
=  \cdots ~0, ~0, ~0, ~x[1], ~0, \cdots$$
creates a response
$$0, x[1]h[0], ~~x[1]h[1], \cdots, ~~x[1]h[n-1], x[1]h[n] \cdots$$
Notice the delay in the response to $x[1]$.  We can continue further
in this vein, but it is best to switch to a more tabular form
and show the various outputs aligned properly in time. We have
$$\begin{array}{l|l|l|l|l|l|l|l}
\text{time} \to & 0 &1 &2 & \cdots & n & n+1 & \cdots \\
\hline
x[0] & x[0]h[0] &x[0]h[1] &x[0]h[2] & \cdots &x[0]h[n] & x[0]h[n+1] & \cdots\\
\hline
x[1] & 0 & x[1]h[0] &x[1]h[1] & \cdots &x[1]h[n-1] & x[1]h[n] & \cdots\\
\hline
x[2] & 0 & 0 &x[2]h[0] & \cdots &x[2]h[n-2] & x[2]h[n-1] & \cdots\\
\hline
\vdots & \vdots & \vdots & \vdots & \ddots & \\
\hline
x[m] & 0 &0 & 0 & \cdots & x[m]h[n-m] & x[m]h[n-m+1] & \cdots \\
\hline
\vdots & \vdots & \vdots & \vdots & \ddots
\end{array}$$
The rows in the above array are precisely the scaled and
delayed versions of the impulse response that add up to
the response $y$ to input signal $x$.

So, we see that if we know the impulse response of an LTI system, we can figure out the response of the LTI system to any arbitrary signal $x$. We don't need any additional information about the LTI system beyond the impulse response.


Z-transform of an impulse signal in discrete time am trying to compute the Z-transform of the following signal
\begin{equation*}
x\left[n\right]\:=\:\sum_{k=-\infty \:}^{\infty \:}\:\delta \:\left[n-k\right]
\end{equation*}
so I thought it would be
\begin{equation*}
X\left(z\right)\:=\:\sum _{n=-\infty \:}^{\infty }\left(\sum_{k=-\infty }^{\infty}\delta \:\left[n-k\right]\right)\:z^{-n}
\end{equation*}
exchanging the order of summation and including the z inside the parenthesis
\begin{equation*}
X\left(z\right)\:=\:\sum _{k=-\infty \:}^{\infty }\left(\sum_{n=-\infty }^{\infty}\delta \:\left[n-k\right]z^{-n}\right)\:
\end{equation*}
then we get
\begin{equation*}
X\left(z\right)\:=\:\sum _{k=-\infty \:}^{\infty }z^{-k}\:
\end{equation*}
which would be
\begin{equation*}
X\left(z\right)\:=\:\frac{1}{\left(1-z^{-1}\right)}
\end{equation*}

Is this correct?? can somebody help me please. Thanks.

There's a ROC (region of convergence) problem in your derivation.

Consider a decomposition of the signal as $x[n] = x_1[n] + x_2[n]$. Then in order for the $Z$ transform of $x[n]$ to exist, the corresponding $Z$-transforms of both $x_1[n]$ and $x_2[n]$ must exist for a set of $z$ which is the intersection of the individual ROCs of the components $x_1[n]$ and $x_2[n]$. If the intersection set is empty then the Z transform does not exist. Applying this to your signal yields.

$$x[n] = \sum_{k=-\infty}^{\infty} \delta[n-k] = u[n] + u[-n-1] = 1 ~~~~,~~~\text{for all n}$$

Where the individual $Z$-transform and their corresponding ROCs are
$$ u[n] \longleftrightarrow \frac{1}{1-z^{-1}} ~~~,~~~ ROC_1 ~~~ |z| > 1 $$
$$ u[-n-1] \longleftrightarrow \frac{-1}{1-z^{-1}} ~~~,~~~ ROC_2 ~~~ |z| < 1 $$

As you can see the intersection of $ROC_1$ and $ROC_2$ is empty; i.e; there is no set of values $z$ for which both of the components has convergent $Z$ transforms. Hence we conclude that the $Z$ transform $X(z)$ of the signal $x[n]=1$ for all $n$ does not exist.


Least-squares digital IIR filter design (with arbitrary responses) I'm studying the IIR filter design that is described in the book: Algorithms for the constrained design of digital filters with arbitrary phase and magnitude responses.

You can get the code at page 171 (at least the main function), and here is an example of a filter design :

M=4;
N=4;
tau=5;
om=pi*[linspace(0,0.2,20),linspace(0.4,1,60)];
D=[exp(-1i*om(1:20)*tau),zeros(1,60)];
W=ones(1,80);
[b,a,e]=mpiir_l2(M,N,om,D,W,0.98);

This is the design of a pass-band linear phase lowpass filter

There is something I don't understand in this code: exp(-1i*om(1:20)*tau) is a way to create regularly spaced points with constant phase change and magnitude = 1. But I don't understand the parameter tau. I tried to change that parameter and I had totally different results: like a notch instead of a lowpass.

How do you choose that number ?
I have to design my filters according to the user desire but how to guess that parameter?
And also how to you choose the phase angles shift ?

Jeff

Well, let's decode it:

om=pi*[linspace(0,0.2,20),linspace(0.4,1,60)];

om is a variable that ranges from 0 to pi with two different spacings (one from 0 to 0.2*pi and one from 0.4*pi to pi.

Then, to find D (the desired response), we do:

-1i*om(1:20)*tau

This takes the first twenty entries in om and multiplies them by tau and -i1 = $-\sqrt{-1}$.

That makes the frequency range (D in the above) from 0 to 0.2*pi*tau (i.e. pi since tau is 5).

Checking page 171, we see:

om frequency grid 0<=om<= pi
D complex desired frequency response on the grid om

which means that om has to be between 0 and pi, and D is the desired response.

If you've changed tau, causing D to have a greater phase (more delay).

What I suspect is happening is that the phase (and, therefore, the delay) you've changed it to cannot easily be approximated by a $4^{\rm th}$ order filter (N or M in your script).


Power spectral density vs Energy spectral density I read the following on Wikipedia:

  Power spectral density:

  The above definition of energy spectral density is most suitable for
  transients, i.e., pulse-like signals, for which the Fourier transforms
  of the signals exist. 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.

I don't quite understand that paragraph. The first part says that "for some signals .. the Fourier transform does not exist".

For which signals (in the context that we are discussing) does the Fourier transform not exist, and we therefore need to resort to the  PSD  rather than using the energy spectral density?
When obtaining the power spectral density, why can't we compute it directly? Why do we need to estimate it?
Finally, on this topic, I have read about methods that use Kayser-windows when computing the PSD over time. What is the purpose of these windows in PSD estimation?

Random process is never ending, non-periodic phenomenon, so taking Fourier transform of its realizations makes no sense, not possible either. However if random process is stationary, then it is for sure that it has some finite power over some band of frequencies. Now, here the question arises that how to compute the power of this stationary random process, (fourier tranform is not possible to be taken directly)? So, what to do? we find the auto-correlation function of the given random process, whose fourier transform always exist. Finally, we take the fourier tranform of this autocorrelation function to get the power spectral density of the given stationary process.

If you integrate the power spectral density of a given stationary process over the interval from -$\infty$ to $\infty$ you ll get the total power contained in the given random process.


Lead compensator vs lag compensator? I already know that lag compensator acts like PI controller and improves steady state and lead compensator acts like PD controller improves transient state but how they achieve their goal?  Despite the fact that apparently both lead and lag compensator have same transfer function?  How are they both different, especially in terms of poles and zeros?

$$ H(s) = {K*\frac{s-z}{s-p}}\\
\\
$$

when |z| < |p|, you have a lead compensator. It adds phase between a certain band, and can help you improve your phase margin. It can also increase your bandwidth. You typically use it to improve your transient as you mentionned. You can think of it as a PD controller cascaded with a low-pass filter. The gain K is optional, but if you set it to $\frac{p}{z}$, the DC gain is 1 which is usefull if you cascade the lead compensator to an existing controller.

I actually never used a PID with a pure derivative in real-life. I have always used a PI controller cascaded with a lead compensator, which is almost the same as a PID but as the advantage of limiting the high-frequency gain, thus improving stability at high frequencies.

When |z| > |p|, you have a lag compensator, which is almost the same as a PI controller. However, you won't have an infinite DC gain. Think of it as a PI controller with a soft reset. I have also seen the term "leaky integrator" used. It could be useful as an alternative to a PI to prevent integrator windup, or if you want a PI controller with a soft reset to slowly reset the integrator.

I have never really used one in real-life. Because, I use appropriate anti-windup strategies. And in the few cases where I needed a "soft reset" for my integrator (for a PLL), I implemented a conditional soft-reset that was only enabled in certain conditions.

If I have time, I will add the Bode plots of both.


Using Principal Component Analysis for corner feature extraction I am on a project of vehicle type Identification (MUV,jeep,Sedan,etc) from images and Here is what I have done up-to now on matlab.

I am working only on the sideview of images.

Input Image :  Masked image

Then I detected the corners easily

Now I got a matrix C with 40 corner coordinates thus I got a matrix of C=[40x2]

I converted it to a linear form using C=C(:); so I got C=[1x80] (as Neural Network destroys coordinates information itself) for 1 image
But I have about 9 image for each class(MUV,Jeep,Sedan) so I made a matrix using C matrix to form a new matrix say MAT=[9x80]

Now I have a corner feature of MAT=[9x80]

1. What Other features could I use or the corner feature would itself be fine enough.(may be SURF features)

2. The Input Format of NN is inMAT=[nxm] where n=number of features and m=total images and target Format is tgMAT=[cxm] where c=number of classes (MUV,Sedan,jeep) and 'm=total images'. Am I right here?

Neural networks can cope with not linearly separable classes, but they're not magic either. They still assume smooth decision boundaries, and I don't think your data would have those.

And before you try to feed the same data to an SVM, RVM, random forest, or any other classifier: I don't think the corner coordinates are good features for any off-the-shelf classification algorithm. The typical input vector in machine learning is something like {age,weight,height,eye color}, where for every input vector, index 2 means "weight". In your input vector, index 2 can sometimes be a corner somewhere near the front wheel of the car, and in another case it's a corner near the roof.

Have you tried simpler shape descriptors, like Fourier descriptors or maybe Hu moments?


$\mathcal Z$-transform, ROC of a system from dividing to others TRUE / FALSE

Given three systems with rational $\mathcal Z$-transform. Systems A and B are not stable with $\mathcal Z$-transform $H(z)$, $G(z)$ respectively. A and B have no common poles. System's C $\mathcal Z$-transform is $F(z)=H(z)/G(z)$, then C is also unstable?

False:

\begin{align}
H(z)&= \frac{1}{1-2z^{-1}}: \textrm{ROC} \quad\lvert z\rvert>2\\
G(z)&= \frac{1}{1-3z^{-1}}: \textrm{ROC} \quad \lvert z\rvert>3\\
F(z)&= \frac{1-3z^{-1}}{1-2z^{-1}}: \textrm{ROC} \quad \lvert z\rvert<2
\end{align}

It is the answer, but doesn't $F(z)$ ROC must be like that of $H(z$), or the division allows to choose the ROC I want?

I agree with your doubts. The $\mathcal{Z}$-transform

$$\hat{G}(z)=1/G(z)=1-3z^{-1}$$

has the ROC $|z|>0$. The ROC of the multiplication of $H(z)$ and $\hat{G}(z)$ must be the intersection of the ROCs of $H(z)$ and $\hat{G}(z)$ (unless there occur pole-zero cancellations, which is not the case here). Since the ROC of $H(z)$ is $|z|>2$, the ROC of $F(z)=H(z)\hat{G}(z)$ must also be $|z|>2$.

Another way to think about it is that since $H(z)$ and $\hat{G}(z)$ are both causal, also $F(z)$ must be causal, i.e., its ROC must be outside a circle in the complex plane.


Update Sub Set of the State Vector in Kalman Filter I question about Kalman Filter:

If I have system state
$$
\mathbf{X} = [x_1\ x_2\ x_3\ x_4\ x_5]^T,
$$
these state elements are independent. I have measurement from a sensor to correct the predicted state. Normally, the correction will perform on all elements of state $\mathbf{X}$ and produce new covariance matrix $\mathbf{P}$.

If I just want to correct the substate only $x_1$ and $x_2$, and do not do anything with $x_3$, $x_4$ and $x_5$, how can I do that ? and how to compute $\mathbf{P}$?

It will depend upon the measurement model. If your sensor provides data for all the states i.e. for $x_1,x_2,x_3,x_4,x_5$, then you can update only $x_1$ and $x_2$ by making the noise covariance corresponding to $x_3,x_4,x_5$ very large. Actually, there will still be some corrections in these three states but they will be small. The higher the covariance corresponding to these three states, the lower is the correction in them.
Alternatively, you can make the last three rows and last three columns of process covariance very small.


DFT/FFT: What did N stand for again? So if I have an image of the size: 1920×1080 and I want to apply a convolution filter of 3x3. Take the formula for DFT $N^2$ and FFT $N \log N$. What should I fill in for N in this case?

Your image is M rows x N cols, and a 2D FFT requires 1D FFTs on all rows followed by 1D column-wise FFTs on the result.

So you have:

$M N \log N + N M \log M = MN (\log M + \log N) = MN \log MN$

operations in total.

The DFT case can be derived using a similar method - this is left as an exercise for the reader.


Explain Shift property of DFT Please check the link

What effect does a delay in the time domain have in the frequency domain?

Here it is said that if you delay your input signal by D samples, then each complex value in the FFT of the signal is multiplied by the constant exp(−j2πkD/N).

My question is that if i have a DFT vector F(u)={11,-3+2j,-1,-3-2j} and D=2,N=4 then I will get {11,3-2j,-1,3+2j}. But if we see this new vector we can see that here no translation
of the signal is done. Signal has not moved to the new location. Just the value of each coefficient has changed in it own location. According to the translation property the answer should be {-1,-3-2j,11,-3+2j}. What is the problem ? Please explain .

But if we see this new vector we can see that here no translation of the signal is done.

You have modified the data in the frequency domain, by multiplying by exp(−j2πkD/N). The translation thus occurs in the time domain. It would be indeed weird if multiplyig a vector by exp(−j2πkD/N) would shift this same vector!

The inverse DFT of [11,3-2j,-1,3+2j] is: [4, 4, 1, 2].

The inverse DFT of [11,-3+2j,-1,-3-2j] is: [1, 2, 4, 4].

The shift property is observed - the signal has been shifted by 2 (modulo 4).


How to recover complex values (downconvert) from FM modulation (composite signal)? How does one recover the complex IQ data from the output of a FM modulator.  I have combined my signals correctly and performed FM modulation on them, now I'm looking for a way to find the IQ values that correspond to each sample output of the FM modulator?

I am writing a modulator for my signal to test the performance of my demodulator and my demodulator depends on complex samples it would receive from a radio.

In other words, how do I perform the down conversion that happens within a tuned radio to produce IQ data?

I am not quite sure if you are asking how to perform the demodulation or how to down-mix a signal at some intermediate frequency to baseband prior to demodulation.  I will assume the prior.

If you have the modulated signal it is quite easy to perform demodulation and recover the original signal.

Simply take the angle (atan2) of the output of a delay discriminator.  If you are concerned with efficiency you can get fancier with your implementation.  Example python code would look like the following:

# Calculate the complex vector between two adjacent data points
vectorDiff = data[1::1] * np.conjugate(data[0:-1:1]);
# Record the angle of the complex difference vectors
fmDemod = np.angle(tmp);

I recently wrote a blog post detailing an FM demodulator I wrote using Python.  Here is a link.  I used this project as an excuse to learn python and numpy.

** EDIT ** Quick and dirty way to perform down-conversion using the hilbert transform prior to FM demodulation (in python)

# Import the plotting library
from matplotlib import pyplot as plot
import scipy.signal
import numpy as np

# Define the setup
fCarrier = 10;
fAudio = 1;
fs = 1000;
timeEnd = 1;
time = np.linspace(0,2,fs*timeEnd);

# Create the signals
carrier = np.sin(2*np.pi*fCarrier*time);
audio = np.sin(2*np.pi*fAudio*time);
audioInt = -np.cos(2*np.pi*fAudio*time);
freqMod = np.sin(2*np.pi*fCarrier*time + 2*np.pi*1*audioInt);

# Downconvert
analyticSignal = scipy.signal.hilbert(freqMod); # wikipedia analytic signal
baseband = analyticSignal * np.exp(-2*np.pi*fCarrier*time*1j); # complex mixing
audioDemod = np.angle( baseband[1::1] * np.conjugate(baseband[0:-1:1]) ); # fm demod


Length of Matlab Raised Cosine Filter Output Complex Vector If I instantiate in Matlab a filter, a vector, and use them in the following way:

vector = (complex(randn(1,k), randn(1,k)))';
filt = comm.RaisedCosineTransmitFilter('RolloffFactor',sdrqpsktx.RolloffFactor,'FilterSpanInSymbols',sdrqpsktx.RaisedCosineFilterSpan,'Gain',1,'OutputSamplesPerSymbol',sdrqpsktx.Interpolation)

  comm.RaisedCosineTransmitFilter with properties:

                     Shape: 'Square root'
             RolloffFactor: 0.5000
       FilterSpanInSymbols: 10
    OutputSamplesPerSymbol: 2
                      Gain: 1
output = filt(vector);

output will be a k*2 x 1 complex double. This is because of the OutputSamplesPerSymbol setting. The documentation on this function states that

  The raised cosine filter has (FilterSpanInSymbols x
  InputSamplesPerSymbol + 1)  taps.

So this filter has 2*10+1 = 21 taps. How then is the output computed? If you take the convolution of the coefficients of filt and vector you will get a (k+taps-1)x1 complex double, which makes sense to me because the output of a FIR filter is the convolution of its impulse response with the input signal. How then is output a k*2 x 1 complex double?

The filter computes OutputSamplesPerSymbol output samples for each input value. You would get the complete convolution if OutputSamplesPerSymbol were equal to $1$ and if you added zeros to your input signal until the filter memory contains no non-zero input values anymore.

This is the standard way a filter routine works; if there is no interpolation or decimation, it computes one output sample per input sample because usually the complete input signal is not stored in a buffer, but it arrives either sample-per-sample or it is split into several blocks.


80-tap Decimation Filter: divide into stages? / Cost of having a number of taps I did a FIR filter hardware implementation with 80 taps (FPGA/VHDL). I will use it later as the decimation filter.
My questions:

Does the number of taps affect a cost of implementation of FIR?

FIR  filters  belong  to  the  class  of  linear  filters,  the  combination  of  N  lower  order  filters  can  create  the desired FIR filter of the higher order. So I was thinking If it makes sense to implement the  combination  of  N  lower  order  filters instead of a one 80 taps filter?

Is the decimation  implemented in one stage usually?

EDIT 1

What is a difference between "filter ( convolution)", " polyphase filter" and "polyphase filter bank"?

Honestly, I dont know ny difference in an implementation.

Does the number of taps affect a cost of implementation of FIR?

Well, insert your definition of "cost", and your question should really answer itself.

FIR filters belong to the class of linear filters, the combination of N lower order filters can create the desired FIR filter of the higher order. So I was thinking If it makes sense to implement the combination of N lower order filters instead of a one 80 taps filter?

Is the decimation implemented in one stage usually?

Indeed, splitting decimation into multiple stages is something that is commonly done, because then only the first filter needs to run at the highest rate.
Now, 80 taps is not what we'd usually call a large filter, so yes, you can do that splitting, it's even advisable.
However, it might really be premature optimization: Your FPGA is probably fast enough to do all 80 MACs at full rate, and it doesn't even have to do that, seeing that you're building a decimating filter, so everything but the most naive implementation would be a polyphase filter implementation, which reduces your complexity to 80/(decimation) per input sample. Seeing that your filter is only 80 taps long, you're probably also not aiming for a very high decimation, anyway.


Does Higher Sampling Rate Increase SNR? I have AM demodulated voice samples that come from a 14bit ADC that I will filter, process, and then send to a CODEC.  I can choose a sample rate of 96kSPS, 64kSPS, 48kSPS, or 32kSPS at the CODEC.  Taking into account most of the voice spectrum is limited to less than 15kHz, would I gain anything by sampling at 96kSPS vs 32kSPS?  I will use either an IIR or an FIR filter between the ADC and CODEC.

For a lower sample rate, the digital filtering would require more taps to acheive the same performance; however, there would be less samples to filter so that trade off seems like a wash.

Sampling at a higher rate will distribute the quantization noise over a wider frequency, thus reducing the noise spectral density due to that quantization noise, with a lot of caveats. For more details on that see
What are advantages of having higher sampling rate of a signal?

In your last paragraph, if you are referring to running the same filter at a lower rate versus a higher rate, at a lower sampling rate the digital filter will require LESS taps for the same filter performance. For more details on that see Filter Order Rule of Thumb and How many taps does an FIR filter need?.

An effective strategy to consider is to use a higher sampling rate which relaxes the analog filtering requirements and can give you more effective bits if needed at the analog to digital boundary (which are realized after subsequent filtering digitally). This is then followed by efficient resampling techniques in the digital domain to get to a lower sampling frequency prior to providing the final filtering as required in your system design (coarse and simpler filters are done at the higher rate, and then final "shaping filters" are done at the lowest rate possible thus minimizing number of taps, power dissipation and resources).

For a further details on decimating to a lower rate for improved filter performance see Fast Integer 8 Hz 2nd Order LP for Microcontroller.


Removing vocal but retaining all background ambient noises I have one audio file containing human speeches, a lot of ambient noises like audience laughter, birds chirping, sounds of natures etc..

Now I want to separate it into two audio files where one should contain all human vocals and other should contain all ambient noises like audience laughter, birds chirping, sounds of natures etc..

I have look into some methods like voicefilter paper from google AI but the result here is the human speech and it's not possible to subtract the resultant with the source sound to get the audio file for other ambient noises, since we can't just assume they are a linear combination.

What you are trying to archive is very hard if at all possible. The pinnacle of research today (including my own) is striving for good results on such tasks. DCASE challanges engage with similar tasks with some interesting results, though the make assumption which cannot be generalized to your case. In DCASE2019 task 3, for example, they assume up to 2 timewise parallel sources where in your case there are no such limitations.

With respect to your conditions and constraints, I would advise looking into source separation and\or classification. If your background audio has a repetitive nature, this paper may help.


Confusion related to norm of dft I have this confusion. I was reading this paper. It says

lets say I have a n dimensional vector X and I take its DFT i.e DFT(X). Then the n Fourier coefficients have the same L2 norm as X. How is it true? I tried on a vector

X = ones(100,1). Applied DFT to X i.e, DFT(X). But norms are not equal

Check the scaling factors used in your DFT/IDFT implementation. There is an arbitrary scaling factor that differs by implementation. Specifically, if you take an input vector $x$ and round-trip it via a DFT followed by an IDFT:

$$
y = IDFT(DFT(x))
$$

Then you may see a couple different results in common implementations:

$y = x$
$y = Nx$, where $N$ is the transform length

Other cases that I wouldn't be too surprised to see:

$y = \sqrt{N}x$
$y = \frac{1}{N}x$
$y = \frac{1}{\sqrt{N}}x$

Which one of these that you see depends on which scaling factors are used in the forward and reverse transforms. If you want the norm of the DFT of $x$ and the norm of $x$ to be the same, then you want the unitary DFT definition:

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


Why does this parallel IIR filter have such a weird form? Referring to "Parallel 1st-order IIR filters" on this wikipedia page, the structure of a 4-parallel IIR filter is discussed. However, the format in which it is implemented seems very odd.

Why is it necessary to iterate the IIR equation N parallel times? Could you not just parallelize it in a direct form. Refer to my chicken-scratch block diagram below.

The main difference between the form the OP is showing and the referenced solution is that the referenced answer isolates the pole and moves it further away from the unit circle resulting in a more stable design applicable to fixed point solutions. In particular the solution the OP gave has delay accumulation from each element that is within the feedback loop (four multiplier and adder elements) which would limit the maximum clock rate that it could operate in a stable condition for in a realizable implementation.

This is depicted in the diagram below showing the synchronous and combinatorial elements. In order for this diagram to accurately match the first order IIR filter (a single feedback tap), the internal functionality shown must be combinatorial. This will run at 1/4 of the clock rate, but will take significantly longer for the state to settle between clock cycles than the implementation shown in the OP's Wikipedia link copied below this image.

OP's Diagram

Wikipedia Diagram
(note this is only showing two of 4 outputs ($y[4k]$,$y[4k+4]$), in Ben's answer he links a pdf that shows the expanded version of this on p. 35 that shows all four outputs that would be needed after each clock cycle. In either diagraml this partial one or the complete one, the feedback loop is the same, which is the main point: it consists of only 1 adder and 1 multiplier running at 1/4 the clock rate, rather than 4 adders and 4 multipliers in the critical feedback loop).


Kalman filter using only location information I have never used a Kalman filter. Is it possible to implement a Kalman filter using only location information?

In addition, if you can implement it, what kind of procedure should be used to program it?
It would be helpful if you could provide a helpful description or site.

If you are new, you might use FilterPy and python. It is such a good start point since it has detailed explanation and manuel.


Can someone let me know linearity, time-invariance, causality, memory characteristic of the system? The input/output system is
$$
\frac{dy(t)}{dt}+2y(t)=2x^2(t)
$$

I want to know this is linear? causal? time-invariant? memoryless?

According to solution, the answer is 'Linear if zero initial conditions, causal, time-invariant, with memory'.

My trial:

$$y(t)=x^2(t)-\frac{1}{2}\frac{dy(t)}{dt}$$

Sorry for my poor drawing.

for input value: $x_1(t)=2x(t)$,
\begin{align}
y_1(t)&=\left(2x(t)\right)^2+\left(-\frac12\right)\frac{d}{dt}\left(y(t)\right)\\
&=4x^2(t)-\frac12\frac{dy(t)}{dt}\\
2y(t)&=2x^2(t)-\frac{dy(t)}{dt}\\
\therefore y_1(t) &\ne 2y(t)
\end{align}

$\because$ not satisfying Homogeneity, Non-Linear.

The system doesn't have $x(t+k)$. (It means $y(t)$ is not made by $x(t+k)$)

$\therefore$ The system is causal because

for input value: $x_1(t)=x(t-t_0)$,
\begin{align}
y_1(t)&=x^2(t-t_0)-\frac{1}{2}\frac{dy(t)}{dt}\\
y(t-t_0)&=x^2(t-t_0)-\frac{1}{2}\frac{dy(t-t_0)}{dt}\\
\therefore y_1(t) &\ne y(t-t_0)
\end{align}

$\therefore$ time-varying.
The system doesn't have any $\displaystyle\int$.

i.e.) $\displaystyle\int x(\tau) d\tau$.

$\therefore$ memoryless

What is the real solution?

There's an error in your 3rd section (linearity) - you replace $y(t)$ with $y_1(t)$ on the left side of the equation but not on the right side. You should have
$$y_1(t)=x^2(t-t_0)-\frac{1}{2}\frac{dy_1(t)}{dt}$$.
Now you see that $y_1(t)=y(t-t_0)$.

As for memoryless - start with your initial equation and integrate both sides. Now you have
$$
\int y(t)dt=\int x^2(t)dt-\frac{1}{2} y(t)dt
$$

and rearranging gives
$$
y(t) = 2 \int x^2(t)dt-2 \int y(t)dt
$$
Does that help?


Modelling zero padded OFDM as circular convolution between channel and FFT frame This is meant to be a puzzle so please preface your answer with spoiler notation by typing the following two characters first ">!"

An often mentioned negative side of the zero padded OFDM is that it does not have the elegant and simple form of a modelled circular convolution between the channel and IFFT frame that a cyclic prefix OFDM provides. However, there is a way that the zero padded OFDM can be modelled as a circular convolution between the channel and signal. How can we acheive this?

Assume an L tap channel and N point IFFT.  $L<<N$.

Also what is the disadvantage to the system, after this modeling can be done.

The answer lies in relationship between Linear Convolution and Circular Convolution. Even Zero padding is not required to model circular convolution using linear convolution. Zero-padding is only required when we are sending OFDM symbols consecutively to cope up with ISI. My point will get clearer in a minute. I will try to show mathematically, by taking an example, where $\vec{d} = d[0], d[1], d[2], d[3], d[4]$ are time-domain samples, meaning IDFT of 5 data symbols in context of OFDM. And the channel is 3 Taps so, $\vec{h} = h[0], h[1], h[2]$. Now, the channel is an LTI system so it takes linear convolution of input $\vec{d}$ with channel impulse response $\vec{h}$, which can be given by the following equation:$$y_{lin}[n] = \sum^{m=2}_{m=0}h[m]d[n-m],$$ignore noise for the sake of this question. Linear convolution output $\vec{y}$ will be of length $5+3-1=7$. Expand the convolution in matrix form to get the following: $$\vec{y_{lin}} = \begin{pmatrix}y_{lin}[0]\\y_{lin}[1]\\y_{lin}[2]\\y_{lin}[3]\\y_{lin}[4]\\y_{lin}[5]\\y_{lin}[6]\end{pmatrix} = \begin{pmatrix}h[0]&0&0&0&0\\h[1]&h[0]&0&0&0\\h[2]&h[1]&h[0]&0&0\\0&h[2]&h[1]&h[0]&0\\0&0&h[2]&h[1]&h[0]\\0&0&0&h[2]&h[1]\\0&0&0&0&h[2]\end{pmatrix}.\begin{pmatrix}d[0]\\d[1]\\d[2]\\d[3]\\d[4]\end{pmatrix}$$This is not circular convolution because the $H$ matrix is NOT Circulant yet. Now remove last 2 rows and add them to top 2 rows to get the following: $$\vec{y_{circ}} = \begin{pmatrix}y_{circ}[0]\\y_{circ}[1]\\y_{lin}[2]\\y_{lin}[3]\\y_{lin}[4]\end{pmatrix} = \begin{pmatrix}h[0]&0&0&h[2]&h[1]\\h[1]&h[0]&0&0&h[2]\\h[2]&h[1]&h[0]&0&0\\0&h[2]&h[1]&h[0]&0\\0&0&h[2]&h[1]&h[0]\end{pmatrix}.\begin{pmatrix}d[0]\\d[1]\\d[2]\\d[3]\\d[4]\end{pmatrix},$$Notice that the last 3 values of $\vec{y_{circ}}$ are unchanged from $\vec{y_{lin}}$, only the first 2 values are changing. And, now the new $H_{circ}$ matrix becomes Circulant. So, the trick lies in removing last $(L-1)$ values of received samples, $y[N], y[N+1], y[N+2], ...,y[N+L-2]$, and adding them to top $(L-1)$ received values, $y[0], y[1], y[2], ..., y[L-2]$, before taking $DFT$ of this new vector $\vec{y}_{circ}$ obtained from original received vector $\vec{y}_{lin}$. Now in context of OFDM, Zero-Padding is required otherwise the last $(L-1)$ samples of this OFDM symbol will get added to the first $(L-1)$ samples of next symbol, and then we will not be able to do the trick which allows us to convert linear convolution into circular convolution. Now, like conventional OFDM RX Chain, DFT can be done on $\vec{y_{circ}}$ to get the following:$$DFT\{ \vec{y_{circ}} \} = H[k].D[k], $$ where $H[k] = DFT\{ \vec{h} \}$. So, basically, the above described $(L-1)$ Zero-Padded system is equivalent to adding Cyclic Prefix in absence of Noise, if we do that extra processing on received samples $y[n]$, in the sense it makes the channel single tap in frequency domain. Disadvantages: First, because we are going to add last $(L-1)$ recieved samples to first $(L-1)$ samples of received vector, hence, in presence of AWGN the noise variance for the first $(L-1)$ data samples will be doubled (Because 2 samples of uncorrelated noise will get added too). However, in absence of noise, both Zero-padding and Cyclic-Prefix will perform same. Second, Zero-padding is anyway required so we are going to waste channel utilisation, but practically speaking, sudden shutdown of TX modulator is not possible. There will be transient and that will ruin the time-domain samples at transmission. And, since TX modulators cannot support complete shutdown of power, zero-padding based OFDM systems are very rarely implemented.


How critical is the selection of the window function in STFTs? I have a sum of periodic signals that I am trying to untangle using time-frequency analysis. I seem to get wildly different results depending on the window length and shape. This is a problem because I want to develop an automated, and hopefully sequential algorithm to do the job.

Window functions have an inherent tradeoff between two of their frequency-domain properties:

Main lobe width: Any tapered window function will cause some "smearing" in the frequency domain. This is visualized by the width of the center lobe in the window function's frequency response. The wider the main lobe, the more difficult it is to resolve two tones that are close in frequency (if they are closer to one another than the main lobe width, they will tend to smear together). So ideally, you would like to have a window function that has a very narrow main lobe.
Maximum sidelobe height: Many window functions have frequency responses that consist of a single main lobe surrounded by repeated sidelobes that decay at some window-specific rate. The height of these sidelobes can make it difficult to resolve two tones that are separated in frequency, but differ greatly in amplitude. So ideally, you would like to have a window function that has very low sidelobes.

The problem: if you decrease the main lobe width of a window function, the sidelobes will grow, and vice versa. So, you need to strike an application-specific balance when choosing a window, based upon the distances in frequency and amplitude that you expect between your signals of interest. Given specific parameters of your system, it's possible to choose a window that (hopefully) meets your requirements.

As far as choosing the length of your window (which is equivalent to choosing the length of the DFT), you're best served with making your observation as long as possible within the constraints that your application might impose (e.g. latency requirements, how long the signals of interest can be considered stationary, computational resources, etc.). Your ability to resolve in frequency is directly proportional to the observation length (measured in time, not necessarily based on the FFT length, which can be zero-padded with no improvement in frequency resolution).


Finding the time domain equation corresponding to a given DFT I'm stuck with the following exercise while self-studying the Discrete Fourier Transform:

  Consider sampling exactly three cycles of a continuous $x(t)$ sinusoid resulting in an 8-point $x(n)$ time sequence whose 8-point DFT is the $X(m)$ shown in figure below. If the sample rate used to obtain $x(n)$  was  4000 Hz,  write the time-domain equation for the discrete $x(n)$  sinusoid in trigonometric form. Show how you arrived at your answer

From $f_{analysis}(m) = {{fs}\cdot{m}\over{N}}   $, I know the "bins" 0 to 7 on the graphics are corresponding to the frequencies in Hz:

$$ \{ 0,500,1000,\boldsymbol{1500},2000,\boldsymbol{2500},3000,3500 \}$$

where the values in bold are corresponding to the non-zero terms.

I also feel like $x(t)$ is of the form:

$$x(t) = A_0\left(cos\left({2\pi}{1500}t\right)+sin\left({2\pi}2500{t}\right)\right)$$

But, (1) I'm not sure my intuition is right (2) I don't know how to prove it (3) I don't see how to calculate the amplitude $A_0$ (4) and finally, is finding $x(t)$ really the answer to the question of writing "the time-domain equation for the discrete () sinusoid in trigonometric form"

(1) I'm not sure my intuition is right

Almost, but you're missing the fact that if you're sampling at 4000kHz, the signals alias; after sampling at 4000kHz, a sine wave at 1500Hz is indistinguishable from one at 2500Hz or 5500Hz, etc.  Also, because the work is being done in complex numbers, negative frequencies have meaning: $e^{-j \omega t}$ is different from $e^{j \omega t}$

  (2) I don't know how to prove it

Go back and study the Fourier Transform more.  You want the Euler identity tattooed onto your frontal lobes: $e^{j \theta} = \cos \theta + j \sin \theta$.  You also want the two offshoots of that: $$\cos \theta = \frac{e^{j \theta} + e^{-j \theta}}{2}$$ and $$\sin \theta = \frac{e^{j \theta} - e^{-j \theta}}{j2}$$

prove these to yourself.  It'll help you out immensely later.

So (after you review the Fourier transform) you have a signal with a real part and an imaginary part: $$\sqrt{32} \left ( e^{j 2 \pi \frac{3}{8}n} + e^{j 2 \pi \frac{5}{8}n} \right) - j\sqrt{32} \left ( e^{j 2 \pi \frac{3}{8}n} - e^{j 2 \pi \frac{5}{8}n} \right)$$  Note that $2 \pi \frac{5}{8} n = 2 \pi (1 - \frac{3}{8}) n$ and recall that $e^{j 2\pi + \theta} = e^{j \theta}$; after some clever trig identities do these look like anything I've just shown you?

  (3) I don't see how to calculate the amplitude A0

It should drop out of the calculations above.

  (4) and finally, is finding x(t) really the answer to the question of writing "the time-domain equation for the discrete () sinusoid in trigonometric form"

Yes, but mind your aliasing!


Non-Linear, Non-Stationary spectral analysis methods! When and where? I have been reading about non-linear non-stationary signal analysis methods and it seems to do this type of analysis the go-to method is the Empirical Mode Decomposition (EMD), then Hilbert Transform (HT) to get instantaneous phase and frequency.

However I have been doing some analysis on a non-linear, non-stationary signal by splitting the time signal into intervals with a Hanning window and Fourier transforming each interval which I think might also be the idea behind the Short Time Fourier Transform (STFT), although I am not sure.

Can anyone tell me what is wrong with this method, and why the EMD / HT method is more advantageous? I have also seen wavelet analysis used in similar cases, is this something I should look into?

Non-linearity and non-stationarity are non-properties. Without more details, they do not say much about the methods that may perform well, and moreover the  choice depends a lot on what you really do: analysis, feature extraction, enhancement, filtering, component separation, restoration? What follows are typical sets of tools you could use:

Your moving-window-Fourier definitely has some STFT flavors, and there are many non-linear extensions: bilinear or quadratic time–frequency distributions (Wigner-Ville, Choi–Williams, Rihaczek)
Non-linearity may come from changing the classical energy preservation into other norms, like robust DFTs (under an $\ell_1$ norm minimization), see for instance Robust $L$-estimation based forms of signal transforms and time-frequency representations, see also Using different algorithms/formulas to create spectrograms
Some methods also alleviate the hypothesis that signal formation is linear: variance-stabilizing transforms,  cepstral analysis or homomorphic signal processing (linear transforms applied in a $\log$ domain)
Generalizations of linear approaches in polynomial or rational fraction modeling, with a recent interest in Rational Optimization for Nonlinear Reconstruction with Approximate $\ell_0$ Penalization.

All of the above can be combined somehow, like EMD at different scales in Intrinsic multi-scale analysis: a multi-variate empirical mode decomposition framework. EMD has many uses and some known limitations at the same time.


Saliency Map based on itti model i want to create a Saliency Map based on itti model via OpenCV(Python). what i'm doing is that:
1- first i change image color space from RGB to LAB color space
2- for each color channel, i create a Pyramid with 6 level
3- from bigger image to smaller, first i make them a fixed size, then differ them via cv2.subtract command.

What i expect is that, in each step of differencing, i should lose details and focus on salient regions but what a get is a black screen. in fact, subtracting removes all image data.

this is my code:

import cv2
import numpy as np

sign =  cv2.imread("sign.jpg")

#size of image
height = np.size(sign, 0)
width = np.size(sign, 1)

#converting LAB color space
lab = cv2.cvtColor(sign, cv2.COLOR_BGR2LAB)

#extracting color channels
l= lab[:,:,0]
a= lab[:,:,1]
b= lab[:,:,2]

a_channel = a.copy()
b_channel = b.copy()

a_prm = [a_channel]
b_prm = [b_channel]

for i in xrange(6):
    a_channel = cv2.pyrDown(a_channel)
    b_channel = cv2.pyrDown(b_channel)

    a_prm.append(cv2.resize(a_channel, dsize=(width,height)))
    b_prm.append(cv2.resize(b_channel, dsize=(width,height)))

a_map = a_prm[0]
for i in range(1,6):
        print(i)
        a_map = cv2.subtract(a_map,a_prm[i])

cv2.imshow('a_map',a_map)

I am developing the same algorithm, but I am using C++ one of the biggest  problem I got was that I  did the math (subtraction and sum) using uchar values obtaining black images. Try transforming uchar to float to do the math and then go back to uchar to visualize the image.


What is the effect on Hilbert transform of signal after multiplication by sinusoid? Question :

Find Hilbert transform of $[u(t-a)-u(t-b)]\cos2\pi f_{0}t\\\\$ such that $\\0<a<b$

my attempt:

we know Hilbert transform of

$[u(t-a)-u(t-b)]\xrightarrow{\mathcal H} \dfrac{1}{\pi}\ln\left|\dfrac{t-a}{t-b}\right|$

but after multiplying left side by sinusoid $\cos 2\pi f_{0}t$  how the right hand side will alter??. Can we expect good results ( like shifting in spectrum by $\pm f_{0}$ which takes place in Fourier transform)

Let $x(t)=\big[u(t-a)-u(t-b)\big]$, and $s(t)=x(t)\cos(2\pi f_0t)$. If we can assume that $f_0\gg 1/(b-a)$, then we can approximate the analytic signal of $s(t)$ (i.e., the signal with all negative frequency components removed) by

$$s_a(t)\approx x(t)e^{j2\pi f_0t}\tag{1}$$

The Hilbert transform of $s(t)$ is the imaginary part of the analytic signal, i.e., we have the following approximation:

$$\begin{align}\mathcal{H}\{s(t)\}&=\textrm{Im}\{s_a(t)\}\\&\approx x(t)\sin(2\pi f_0t)\\&=\big[u(t-a)-u(t-b)\big]\sin(2\pi f_0t)\tag{2}\end{align}$$

The same result is easily obtained by applying Bedrosian's theorem, which assumes that $x(t)$ and $\cos(2\pi f_0t)$ have no spectral overlap. This is approximately satisfied if the same condition as above is assumed to be true: $f_0\gg 1/(b-a)$.

This condition basically says that the signal $s(t)$ is a band pass signal, which implies that $x(t)$ is a low pass signal. The Hilbert transform of such amplitude modulated signals always equals the original low pass signal $x(t)$ multiplied by the Hilbert transform of the carrier:

$$\begin{align}\mathcal{H}\big\{x(t)\cos(\omega_0t+\phi)\big\}&=x(t)\mathcal{H}\big\{\cos(\omega_0t+\phi)\big\}\\&=x(t)\sin(\omega_0t+\phi)\tag{3}\end{align}$$

The exact Hilbert transform of the given signal can only be computed by evaluating the integral

$$\begin{align}\mathcal{H}\big\{s(t)\big\}&=p.v.\frac{1}{\pi}\int_{-\infty}^{\infty}\frac{s(\tau)}{t-\tau}d\tau\\&=p.v.\frac{1}{\pi}\int_{a}^{b}\frac{\cos(2\pi f_0\tau)}{t-\tau}d\tau\tag{4}\end{align}$$

but I don't think that that was the idea of the exercise.


What does zero-mean noise mean? I don't know where to start from. and this question might be very silly,
I have found such as comment from someone at Poisson noise and curve fitting - denoise first?,
"I think if your noise is zero-mean it doesn't matter, if not it does (you should de-noise first). I'm not entirely sure though so I won't put it in an answer yet."
Q1. Actually, I want to know what does zero-mean noise mean? and what does it not matter means?
Q2. Why we presume that the noise in Gaussian filter have sigma 1 and mean 0 ?

It is an easy concept but a sometimes it is hard to express it compactly.
The difficulty arises because of additional concepts to be known for any formal description.

As you may know, mathematical analysis of noise variables and signals are handled via the tools of probability theory, in particular by the stochastic (random) processes.

The usual progress is such that first properties of random variables are described and then properties of random signals are derived based on it.

The mean for a random variable, X, is defined to be its Expected Value E(X) and computed based on its Probability Density Function $f_X(x)$ such as $\mu_X = \int {xf_X(x) dx}$ with the limits from minus to plus infinity for a continuous random variable.

So a zero mean random variable is that one for which the above integral is zero. (for a discrete random variable the integral is replaced by a sum and pdf is replaced by a pmf). This is also called as the average value.

Random signals are considered to be made up of an infinite set of random variables each making a single value at a single instance of observation of the instance being created. Hence there is the concept of ensemble and time averages.

Now coming to random signals X(t) (noise) expected value of a random signal is also expressed as E(X(t)) which for a Stationary (or at east Weakly Stationary up to second order) Process is a fixed value.

Your answer is that a zero mean noise is that one for which E(X(t)) is zero for all t.

However above is a theoretical description of Mean. In practice to compute mean and other theoretical parameters, you will refer to statistical computation which deal with practically observed data and computations based on that.

The key which connects statistically made (practical) computations to theoretical computations is the Ergodicity Theorem which states that for an ergodic random process theoretical computations are equal to statistical ones.


Sampling frequency and oversampling factor A quite naïve question. I need some clarification regarding choosing the sampling frequency and oversampling factor. The scenario is as follows:

I have an OQPSK modulated sequence with symbol rate 2 M Symbols/sec. In order to transmit this through an AWGN channel, I am trying to half sine pulse shape this modulated sequence. I have read from a couple of online resources that before pulse shaping the sequence should be oversampled. I am currently confused by the factor that I need to oversample the sequence.

Regarding the sampling frequency, at Baseband I am confused as to what I should take as the sampling frequency. I assume that at passband of 2.4 GHz, I am considering the sampling frequency of 6 GHz but I am not sure if this is right and how should I decide on the sampling frequency if I do not consider the passband and just the baseband frequency.

The rule you must remember when deciding on your sampling frequency is: it must be at least double the highest frequency signal that will ever appear in your simulation.

Consider your baseband signal. The highest frequency signal is the half-sine pulse (its bandwidth is infinite in theory). We know, from the Fourier series of a full-wave rectified sine, that its spectrum decays with frequency as $$\frac{1}{\pi(4n^2-1)}.$$ I would go with $n=3$, which corresponds to a frequency of $3f_0$. If your pulse rate is 2 Mpulses per second, then $f_0=2$ MHz, the pulse bandwidth is 6 MHz, and your sampling rate must be 12 MHz.

Regarding passband: If you must simulate this system at a high carrier frequency, you know that the bandwidth now will be 2.4 GHz plus 6 MHz, or 2.406 GHz, so your sampling frequency needs to be double that. However, before doing this, please learn about the complex envelope. There is nothing to be gained by simulating the passband system, and you can lose a lot of time.


The magnitude spectrum of a sharpening filter I'm trying to derive an expression for the magnitude spectrum of the following sharpening filter.
$$
g(m,n) = \delta(m,n)+\lambda (\delta(m,n) - h(m,n))
$$
where $\lambda$ is some positive constant and $h(m,n)$ is a low pass filter.
$$h(m,n) = \begin{cases}
\frac{1}{25}, & \text{for $|m| \leq 2, |n| \leq 2$} \\
 0, & \text{otherwise}
\end{cases}$$
Using the linearity of Fourier transform, I get
$$ G(u,v) = 1+\lambda(1-H(u,v))$$
where
$$
H(u,v) = \frac{1}{25} \cdot \frac{\sin\frac{5\mu}{2}\sin\frac{5\nu}{2}}{\sin\frac{\mu}{2}\sin\frac{\nu}{2}} \cdot e^{j\cdot 6(\mu+\nu)}
$$
So, to get the magnitude spectrum, I make the following steps
\begin{align*}
|G(u,v)|^2 &= \big[(1+\lambda)-\lambda H(u,v)\big] \big[(1+\lambda)-\lambda H^*(u,v)\big] \\
& = (1+\lambda)^2-(1+\lambda)\lambda\cdot 2\Re \big\{H(u,v)\big\}+\lambda^2|H(u,v)|^2
\end{align*}
where
$$
\Re\big\{H(u,v)\big\} = \frac{1}{25} \cdot \frac{\sin\frac{5\mu}{2}\sin\frac{5\nu}{2}}{\sin\frac{\mu}{2}\sin\frac{\nu}{2}} \cdot \cos(6(\mu+\nu))
$$
Therefore, to obtain the magnitude spectrum I should take the square root of the expression above.
However, the picture I get when I'm trying to plot it does not coincide with something that I'm expecting to see. I mean it's a high pass filter, hence there should be something like that.

But I get this.

Maybe there is a way to simplify this expression or maybe there is something that I'm doing wrong.
Thank your for your help!

UPD: I've figured it out myself. I was inattentive with the exponentials while deriving the expression for the low-pass filter $H(u,v)$. There shouldn't be $\cos(6(\mu+\nu))$ term, then the result is fine and $$\Re \big\{H(u,v)\big\}=\frac{1}{25} \cdot \frac{\sin\frac{5\mu}{2}\sin\frac{5\nu}{2}}{\sin\frac{\mu}{2}\sin\frac{\nu}{2}}$$
Therefore,
\begin{align*}
|G(u,v)|^2 &= \big[(1+\lambda)-\lambda H(u,v)\big] \big[(1+\lambda)-\lambda H^*(u,v)\big] \\
& = (1+\lambda)^2-(1+\lambda)\lambda\cdot 2\Re \big\{H(u,v)\big\}+\lambda^2|H(u,v)|^2 = [1+\lambda - \lambda H(u,v)]^2
\end{align*}
And the result is exactly what I was expecting to get.


What audio effects, filters, distortions, etc. create a "vinyl" effect? I want to make any song sound like it's being played on a gramophone record in the 1920s. What specific filters would I apply to make this happen?

I'm looking for technical details, not some magical program that does it for me. I'm trying to create a programming interface that does this, so a user interface doesn't work. I'm a Math major with understanding of FFTs, but I don't know all the lingo.

I don't know the precise formula for 1920's gramophone, but this is the general process for audio antiquing and should get you started. Depending on the exact settings, you can get vinyl, telephone, etc.

Convert to mono.
Add white noise.
Bandpass. I would start with a HP at 80 HZ and a steep LP at 8 kHz, but that's a guess. You'll have to use your ears until it sounds like the gramophone. You should try Chebychev and Butterworth filters if available.

You may also want to add "clicks" and "Pops". You can model those in a variety of ways. One way to start is with a Dirac function and filter it with bandpass functions. Depending on how strict your bandpass in step three is, you might be able to just add the clicks and pops as Diracs at the same time you add white noise, but then all your clicks and pops will sound about the same, which isn't very realistic.

Some other things you may want to try for more realism:

applying compression and saturation after step 1.
applying wow and flutter emulation either before step 2 or after (or, for most accuracy, two stages of white noise before and after the wow and flutter emulation).

You may want to ask on https://video.stackexchange.com/ to see exactly what filter settings they would recommend for 1920's gramophone in particular.

For tips on implementing basic audio filters, I have a blog post: http://blog.bjornroche.com/2012/08/basic-audio-eqs.html


What happens with signal in frequency spectrum when it is time shifted in time spectrum? I have some trouble to understand what is going on with signal in frequency spectrum when it is time shifted in time spectrum.

I am hoping that somebody will help me to understand that.

Thanks you very much.

Each frequency in the FT of a time shifted waveform is rotated in phase by an amount proportional to the frequency and proportional to the amount of time shift.

If you delay a pure sinusoid by 25% of its period, it's phase referenced to any fixed point in time will change by pi/2 radians.  Delay a slightly higher frequency sinusoid by the same amount of absolute time, and it's phase will change more.  So the phase change will change with frequency.

If looking at a 3d plot of the FT of a time shifted signal, it looks like taking the FT before the time shift and twisting it.  The twist will be linear, e.g. proportional to frequency.  The more time shift, the greater the amount twist per unit of frequency (or twist revolutions per graph width).


How do I get the most accurate camera calibration? First off, I hope this is the correct Stack Exchange board. My apologies if it is not.

I am working on something that requires me to calibrate the camera. I have successfully implemented the code to do this in OpenCV (C++). I am using the inbuilt chessboard functions and a chessboard I have printed off.

There are many tutorials on the internet which state to give more than one view of the chessboard and extract the corners from each frame.

Is there an optimum set of views to give to the function to get the most accurate camera calibration? What affects the accuracy of the calibration?

For instance, if I give it 5 images of the same view without moving anything it gives some straight results when I try and undistort the webcam feed.

FYI to anyone visiting: I've recently found out you can get must better camera calibration by using a grid of asymmetric circles and the respective OpenCV function.

You have to take images for calibration from different points of view and angles, with as big difference between angles as possible (all three Euler angles should vary), but so that pattern diameter was still fitting to camera field of view. The more views are you using the better calibration will be. That is needed because during the calibration you detect focal length and  distortion parameters, so to get them by least square method different angles are needed. If you arn't moving camera at all you are not getting new information and calibration is useless.
Be aware, that you usually need only focal length, distortion parameters are usually negligible even for consumer cameras, web cameras and cell phone cameras. If you already know focal length from the camera specification you may not even need calibration.

Distortion coefficient are more present in "special" cameras like wide-angle or 360°.

Here is the Wikipedia entry about calibration.
And here is non-linear distortion, which is negligible for most cameras.


How do you sub-sample an FFT? I'm currently trying to understand this paper on a Sparse Fourier Transform.

On page four, there is a section on a sub-sampled FFT. The purpose of this section is to show that you can compute a sub-sampled FFT using less computations than the original FFT. So if I have a signal of size N, I could compute a subsampled FFT of size B which would require O(N + B log B) operations.

However, I'm confused as to how you would obtain this subsampled FFT. The idea is that we want the original FFT spectrum to be sampled at B locations. Obviously we could miss coefficients by only sampling the spectrum at B locations, however, it seems like the bigger problem is that subsampling the signal would result in aliasing frequencies.

Essentially I'm not sold on the fact that the described algorithm would result in a sub sampled FFT of the original spectrum. Instead, I think the spectrum could be aliased resulting in incorrect coefficients.

it seems like the bigger problem is that subsampling the signal would result in aliasing frequencies

No, subsampling in frequency domain corresponds to aliasing in time domain. So the idea here is to purposefully alias in time domain so that you get a sub-sampled FFT. That is exactly why 'Claim 3.7' of paper mentions $y_i=\sum_{j=0}^{n/B-1}x_{i+Bj}$. These are $n/B$ copies of $x[n]$ shifted by $B$ and overlapped. Each copy is positioned at $0,B,2B,\ldots,(\frac{n}{B}-1)B$. This aliased signal can be used to compute FFT ($B\log_2 B$ operations) to arrive at subsampled FFT. Total number of operations will the operations need to create the aliasing (summing of non-zero values of $x$ which equals $O(\text{supp}(x))$ in addition to the operation needed to compute FFT of aliased $x$ ($O(B\log_2 B)$).

For example, if a time domain signal is $1024$ samples long, ideally you need $1024$ point FFT. But if you sub sample by taking every 8th sample, $B=128$, $n/B=1024/128=8$, the time domain signal at $\tilde{x}[0]$ will be having addition from $x[0]$,$x[128]$,$x[256]$,$\ldots$,$x[896]$. Like this for $\tilde{x}[n]=\sum_0^{7}x[128k+n]$. You can now use this $\tilde{x}$ to compute $B$-point FFT to arrive at the Subsampled FFT.


Upsampling before performing pulse shaping I know what an upsampler does. We insert L-1 zeros between the samples. As a result, it creates spectral repetitions. What if we have a symbol impulse stream and then we apply a Root Raised Cosine Filter, this would automatically perform pulse shaping. Why do we need to upsample that symbol stream to cause interpolation? The Pulse shaping filter is doing that.

Please explain intuitively if possible because I have been stuck on this topic for a while. Thank you in advance!

The reason is that you want to space the raised-cosine pulses by the symbol interval, $T_p$. Consider the signal you want to create: $$s(t) = \sum_k a_k p(t-kT_p),$$ where $p(t)$ is  your prototype raised-cosine pulse, and $a_k$ are the symbols. Notice that pulses are spaced $T_p$ seconds.

To re-create this with a filter, you can write the equation for the signal as follows: $$s(t)=\sum_k p(t) \ast a_k\delta(t-kT_p),$$ where $p(t)$ is interpreted as the filter's impulse response. Notice that the impulses are spaced $T_p$ seconds, and this results in the raised-cosine pulses being spaced by the same amount.

If you want to simulate this process in discrete-time in Matlab or similar tools, then you can create a vector p for $p(t)$, and a vector d for $a_k\delta(t-kT_p)$. Let's say your symbols are stored in vector a. Then, you can create d by upsampling a in such a way that there are $T_p$ seconds between symbols:

d = [ a[1], 0, ..., 0, a[2], 0, ..., 0, a[3], 0, .... ]

You need to insert as many 0s as necessary, according to your sampling frequency, to have $T_p$ seconds between symbols.

Then, you can generate the vector s corresponding to $s(t)$ with s=conv(p,d).


window for short time fourier transform When we implement the 1D-short time fourier transform, the formula is as such:

$$
Sf[m,l]=\sum_{n=0}^{N-1}f[n]g[n-m]\exp\left(\frac{-i2\pi l n}{N}\right).
$$

the window $g$ used is normally symmetric, why is this so? (Hint: Use the phase of Short time fourier transform.

You are talking about symmetry in two contexts :-

Gaussian
Window

I will talk about these two contexts one by one

Gaussian

A symmetric 2-D gaussian (I assume 1-D as in 1-D the question of symmetry does not arise) means that it does not have any directional preference (dictated by equal variances $\sigma_x$ and $\sigma_y$) . Now the job of Gabor Transform is to give us a joint time-frequency distribution (i.e telling us about frequency and phase content of a local image patch)
Now consider the following two cases :-

Symmetric Guassian ($\sigma_x$ = $\sigma_y$)
In this case the gaussian does not have any directional preference.

In the above image, suppose you are centered on the central pixel and consider a patch around it which is indicated in the above image by a black rectangle in the center. Following points are important for this image :

Doing a Gabor Transform over that patch would yield a time-frequency representation, which involves all the pixels in the patch accordingly weighted by the symmetric gaussian.
All the directions at a certain distance from the origin get the same weigthage.

However, if a gaussian oriented at an angle of $135^{\circ}$, then :-

the white principle diagonal pixels would get almost all the weightage while the other pixels will get little weightage.
The time-frequency representation will be mainly a time-frequency representation around the principle diagonal.
If the values of $\sigma_x$ and $\sigma_y$ are sufficiently small, the gabor transform will be very very close to the time-frequency representation for the straight line which marks the principle diagonal.

Thus, for the Gaussian function,

Symmetric gaussian is non-directional.
Non-Symmetric gaussian is directionally oriented.
Non-Symmetric gaussian will be useful in cases where the time-frequency analysis of a directionally oriented feature in a patch has to be carried out.

Window

The Window size basically represents the time-resolution (or space-resolution in the case of 2-D signals like images). A narrow window means that the focus on some local patch of the image is high and hence a more faithful time resolution of the Time-frequency analysis will be possible.

However, this is limited by the Heisenberg's Uncertainty principle. Although, Gabor Transform attains the theoretical lower bound for the uncertainty principle, it stays within the bounds of the principle. Thus, changing the window size is akin to changing your time resolution.

A Non-symmetrical window (like $3\times7$) is generally not chosen because generally we do not have a reason or a preference for changing our space-resolution in one direction (x or y) relative to other. In case there is a concern for the same, it can be done. However, I have never till now seen such a situation.


Unknown symbol/expression in text about adaptive filters (cst) I am currently reading a chapter about adaptive filters from the Springer Handbook of Speech Processing.

In a formulation of the variable stepsize normalized least mean squares (VSS-NLMS)-algorithm, i found an expression that I do not know yet, and searching the internet did so far not yield results i can make sense of.

The algorithm is expressed as

$$\mathbf{\hat{h}}(k) = \mathbf{\hat{h}}(k-1) + \mu_\text{NPVSS}(k)\mathbf{s}(k)e(k),$$

where $\mathbf{\hat{h}}(k)$ is the current estimate of the desired filter impulse response, $\mathbf{s}(k)$ the source signal and $e(k)=x(k)-\hat{x}(k)$ the error difference between filter output and actual output (using the real $\mathbf{h}(k)$).
The stepsize for this algorithm is defined as
$$\mu_\text{NPVSS}(k) = \frac{1}{\mathbf{s}^\text{T}(k)\mathbf{s}(k)} \left[ 1-\frac{\sigma_b}{\sigma_e(k)} \right],$$

with $\sigma_b$ and $\sigma_e$ denoting the square root of the variance/power of a noise process $b(k)$ uncorrelated to the source signal and of the error difference, respectively.

For a practical implementation it is suggested to compute the fraction $[\mathbf{s}^\text{T}(k)\mathbf{s}(k)]^{-1}$ as $[\delta + \mathbf{s}^\text{T}(k)\mathbf{s}(k)]^{-1}$ instead.

In this context, the regularization parameter is defined as
$$\delta = \text{cst} \cdot \sigma_s^2.$$

I am not too familiar with the concept of regularization, so I have no clue what the term $\text{cst}$ means in this algorithmic formulation.

What does it mean? Any pointers or literature recommendations would be very much appreciated.

The standard normalized step-size LMS algorithm computes the current step-size according to

$$ \mu = \frac{c}{s_k^T \cdot s_k}  $$

where $c$ is a suitable scale factor and $s_k^T \cdot s_k$ is the total energy of the current tap inputs. The algorithm aims to adjust step size according to input signal power; when input has large power then decrease the step-size so that filter coefficient update will not deivate much. But if the signal goes to zero or very small levels (silence durations, or small variance sections), then the division will possibly be very large and algorithm may produce unbounded (very large) outputs or even diverge...

To prevent this from happening, a small nonzero factor; the regularisation factor, is added to the denominator

$$ \mu = \frac{c}{s_k^T \cdot s_k + \delta} $$

to limit the maximum of the fraction to $c/\delta$ when input power goes to zero.


What is the RFFT and iRFFT complexity? What is the RFFT and iRFFT complexity for 1D and 2D?
Here is some starting point.
1D RFFT

The FFT is $\mathcal{O}(N\log N)$
RFFT computes $N//2 + 1$ frequencies instead of $N$.

Thus, RFFT is $\mathcal{O}([N//2 + 1]\log N)$ = $\mathcal{O}\left (\frac{N}{2}\log N \right )$

Question: What about iRFFT? Is it $\mathcal{O}(N\log N)$ or $\mathcal{O}\left (\frac{N}{2}\log N \right )$ ?

2D RFFT

The 2D-FFT is $\mathcal{O}(2N^2\log N)$
RFFT computes $(N//2 + 1)N$ frequencies instead of $N^2$.

Thus, RFFT is $\mathcal{O}([N//2 + 1]N\log N + N^2 \log N)$ = $\mathcal{O}\left (\left [ \frac{3}{2} N + 1\right ] N \log N \right )$ = $\mathcal{O}\left (\frac{3}{2} N^2 \log N \right )$

Question: What about iRFFT? Is it $\mathcal{O}(2N^2\log N)$ or $\mathcal{O}\left (\frac{3}{2} N^2 \log N \right )$?

Am I getting anything wrong?

In Big-O notation, $\mathcal{O}(N\log N)= \mathcal{O}(\frac N2\log N)$, as big-O notation only describes an upper bound for growth up to a constant factor, so your first question is kind of impossible to decide; in practice, the effort for a cleverly implemented real-valued FFT will be less than that of a (normal, complex) FFT.
But whether it will be exactly half might be doubted: as basically everywhere in numerical algorithm, the number of operations isn't actually sufficient to describe the effort or time needed. When I look at my PC, it does about 200 to 1000 multiplications in the time it needs to access a random address in RAM that's not cached. So, you'll see big jumps in actual time when your data size crosses cache size boundaries, which might happen later with real-valued data; or it might not matter at all.
The same problem, big-O notation only requiring the existince of a constant, finite factor $k$ such that the complexity of the algorithm stays below $k$ times the function in the parentheses, remains for your 2D question.
Generally, you'll find that the pure description through numbers of multiplications and additions makese nearly no sense, unless you're dealing with very large numbers, so that you can be asymptotically sure to never be able to e.g. align memory accesses. I wish people would teach  that more intensely when they tell people about algorithmic complexities: for algorithms where the input data is almost certainly not randomly spread across memory, the memory access usually becomes a very dominant factor!
In the case of FFTs, the problem is even more prominent: Everyone (and their grandma) knows that the $N$-point FFT's complexity $c(N)$ is in $\mathcal O(N\log N)$, but nobody notices that there's quite possibly a very large constant factor in $\lvert c(N)\rvert \le k\cdot N\log N$ (and that the base of the logarithm isn't even defined), and that the formal definition especially notes that this only needs to hold true for all $N$ larger than some very large $N_0$.
But it gets worse: Practically, we know that on "normal" computers, on most FPGA IP cores and DSP slices, and most other platforms, power-of-2 FFTs are very fast, whereas FFTs with larger prime factors do still exist. But, I wouldn't be surprised at all if the $N_1=4\cdot 2^{n}$-FFT was faster than the $N_2=3\cdot 2^{n}$-FFT, even if by judging the O-notation, for sufficiently large $n$ that shouldn't be the case! People really mislead when they tell other people the FFT is always in $\mathcal(N\log N)$; you'll need to put the FFTs with the same prime factors in separate classes and only change the exponent of one of these prime factors for that to actually be a sensible statement.
Sorry that I'm such a downer about big-O notation – it's just that beyond showing how a Cooley-Tukey / Radix-2 FFT can save you computation for power-of-2 FFTs, it's very rarely used to imply anything correct, less frequently even to arrive at a useful result!


How should we read FIR filter phase response graph? What should we see from the FIR phase response graph, as an example for a linear phase? And what does this straight line at $0.5\pi$ mean?

First of all you see that the phase is a piecewise linear function, so it's a linear phase FIR filter. There's a phase jump at half the Nyquist frequency, which shows that the filter has a zero at that frequency. Note that the phase jumps by $\pi$ corresponding to sign inversion.

You can also see what type of linear phase FIR filter it is. Since the phase is zero at $\omega=0$ it must be either type I or type II (i.e., even symmetry). If it were type III or type IV (odd symmetry), the phase would jump at $\omega=0$, and the phase at $\omega=0^+$ would equal either $\pi/2$ or $-\pi/2$. You can distinguish type I and type II by looking at the phase at the Nyquist frequency. A type I filter (odd filter length) must have a phase of either zero or $\pi$ at Nyquist, whereas the phase of a type II filter jumps at Nyquist because a type II filter always has a zero at Nyquist. You can't see the phase jump here (because the figure doesn't go beyond Nyquist) but you can see that the phase equals $-\pi/2$ just before Nyquist. A type II filter's phase always equals $\pm\pi/2$ at $\omega=\pi^-$, and it jumps by $\pi$ at $\omega=\pi$.