Trying to apply Parks-McClellan optimal FIR filter design My goal is to apply the Remez exchange algorithm in R to design an FIR band pass filter for a specific frequency window. The question is not supposed to be referring to R of course but about how to choose the parameters correctly and why the result looks like as it does.

The toy signal is an additive combinations of four sine waves each adding a single frequency component to the final signal. The frequencies are 1/20, 1/300, 1/5000 and 1/7000. Its FFT plot is shown in figure (s).

Now what I apparently didn't get yet is how I would have to choose the parameters for the filter design - being the pass band and the stop band cutoff points.

The 1st filter uses the cutoff points:

0
1/10 - 2/1000
1/10 - 1/1000
1/10 + 1/1000
1/10 + 2/1000
1

With 1 being the Nyquist-frequency and the valid cutoff point interval being (0,1).

As you can see in figure (r1) this results in an attenuation of all frequency components except for 1/5000. (I would assume that also the frequency component 1/5000 is actually attenuated but relatively less then the other components)

So my first set of questions would be - why is that so? How could I calculate f.x. the frequency located in the center of the pass band? What is the quantitative relation between cutoff points on an interval (0,1) and the pass band?

The 2nd filter uses the cutoff points:

0
1/10 - 2/500
1/10 - 1/500
1/10 + 1/500
1/10 + 2/500
1

And the result is an (relative) attenuation of all components except for 1/300. What I am wondering about in this case is why the pass band is apparently completely different from the above filter even though I just made the filter's pass band wider. I would assume at least that 1/5000 is still located within the pass band and not attenuated. A possible explanation might be that 1/300 is now included within the pass band and relative to 1/5000 less attenuated making 1/5000 look just smaller.

I would like to add that I am - as you can see - quite new to this subject and actually from a different area and trying my best to understand what's going on.

The canonical question would be: How do I have to choose the cutoff points in case of a possible interval (0,1=Nyquist-Freq) to attenuate all frequency components except for 1/N and its close neighbourhood?

library(signal)
par(mfrow=c(5,1))

# the toy signal
ch <- sin(2*pi*1:100000/floor(100000/20)) +
  sin(2*pi*1:100000/floor(100000/300)) +
  sin(2*pi*1:100000/floor(100000/5000)) +
  sin(2*pi*1:100000/floor(100000/7000))

# (s)
barplot(abs(fft(ch))[1:10000])

filter <- function(c0,d1,d2) {

  # the filter design with cutoffs specified
  fir <- remez(n=1000,
    f=c(0,c0-d2,c0-d1,c0+d1,c0+d2,1),
    a=c(0,0,1,1,0,0))

  # (fN)
  plot(freq$f,abs(freq$h),type="l")
  freq <- freqz(fir,n=100000)
  y <- signal::filter(as.vector(fir), 1, x=ch)

  # (sN)
  barplot(abs(fft(y))[1:10000])
}

filter(1/10,1/1000,2/1000)
filter(1/10,1/500,2/500)

There are a few problems with your designs. First, the desired bandwidths and transition bands (between passbands and stopbands) are very narrow. Second, some of the center frequencies are very low, especially the ones at 20/100000 and 300/100000. And finally, your chosen filter length is very large. I do not know the implementation of the Remez exchange algorithm that you use, but I've seen several implementations that run into serious numerical problems for such filter lengths. You can see in your designs that around Nyquist there are some small peaks in the filter response. These are due to numerical problems.

So you should try to make the passband and the transition bands a bit wider and reduce the filter length to see if you can get a reasonable design. From there you could try to make the specifications a bit more stringent and increase the filter order and see how far you can go.

In Octave, the following specification gives a reasonable result (center frequency: 0.1, filter length: 300, '1' corresponds to Nyquist):

h=remez(300,[0,.07,.085,.115,.13,1],[0,0,1,1,0,0]);


What is the physical interpretation of Lowpass/Highpass filtering? In Image/Signal processing context we have seen/we know that if there are say 4 samples and if we take an average of those 4 samples, then we say that the result sample is a low pass filtered output sample. This is very relevant in Image processing context. Now then, what is the operation which implies a high pass filtering, is it finding some max of the 4 samples and putting that as the output sample, or what other mathematical operations implies high pass filtering as averaging implies low pass filtering.

The example you gave of taking 4 samples and taking the average of it is sort of a poor-man's low-pass filter. Generally things aren't as simple as that. But for understanding sake there is some value in using these simple examples.

A low pass filter is indeed like taking 4 samples and taking an average of it. Ex:

samples = [6 1 -10 -4];
avg_value = mean(samples) = -1.75

The high pass filter is removing the "DC" term. Or more generally, it is removing the data that isn't changing. A simple way of thinking of this is to subtract your avg_value from every sample. Ex:

samples = [6 1 -10 -4];
avg_value = mean(samples) = -1.75;
high_pass = samples-avg_value;
high_pass: [7.75 2.75 -8.25 -2.25]

Now if you take the average of the "high passed" signal you end up getting 0.

These two 'filters' that I mention are both to the extreme in that you have one filter that only gives you DC and another filter that only removes DC. Basically what you end up getting is this ideal filter where the low-pass filter gives you the green and the high pass filter gives you the yellow.

Most of the filters you will use will have a response that looks more like this for a low pass:

and this for a high pass:


Filter design with 8-tap latency on recursion I'm curious about the feasibility of designing a noise shaping filter (but might generalise to any recursive filter) with the constraint that the most recent output samples aren't available for several iterations of the filter.

The use case I have in mind is reducing the word length of an audio stream sampled at 48kHz, and the filter would try to keep most of the quantisation noise out of the 200Hz-8kHz band (where human hearing is more sensitive); but with the additional constraint that the most recent seven outputs are not available to the filter.

I'm applying this constraint so that the output can be produced several samples at a time using a mixture of SIMD and instruction-level parallelism -- working with a notional four-lane SIMD and two-fold unrolling.

There's a lot of work for each sample which could be performed in parallel (scaling, quantisation, dither) if not for the recursive factor.  Normally one might hope to use SIMD for the multiply-accumulate of the coefficients, but that's not a great proportion of the work if the filter has only a handful of taps.

So I'm wondering if it's actually possible to design a filter that doesn't get in the way, here, but is still effective for its intended purpose.

Stereo would obviously offer an opportunity to halve the delay, but avoiding that assumption is preferable, as it makes mono a failure case and multi-channel becomes complicated.

This is tricky. Any filter with 8 tab latency between coefficients can be represented in the z-domain as a rational function in $z^{-8}$ . That's basically the same as saying the impulse response has 7 zeros between each non zero tap. Any filter with this property is periodic in the frequency domain. You can pick any transfer function from -pi/8 to pi/8 but the regions from pi/8 to 3*pi/8,etc. will just be repetitions of that.

You may be able to do better with varying delays but this is a complicated design problem


Interference Distribution Assumption in NOMA IoT Scenario I am assuming a NOMA IoT Uplink scenario, which is interference marred. I have to ascertain the SINR at the Access Point (AP) for the NOMA users. A sub-channel is allocated to various users and interference is considered at each sub-channel. What is the realistic and best assumption for Interference Distribution at each sub-channel?

The interference intensity is different at each sub-channel but varies with same distribution.

interference intensity is different at each sub-channel but varies with same distribution

You're already using the term intensity, so this sounds like you're most likely implicitly using a Poisson Process, which is a stochastic typical model for "arrival processes", e.g. packets appearing in a network. It's based on the idea that the probability of an event occurring within a time period is simply proportional to the length of that period, and does not depend on the start time of the period.

Really check your literature thrice for them indications of where the term "intensity" came from. It's totally possible that something in there confirms the above (i.e. they say that the intensity $\lambda$ defines the factor between expected number $N$ of occurences and the duration of observation $T$, $\mathbb E[N(T)] = \lambda T$); it's also totally possible that other assumptions make this model impossible ("we assume interference to happen in bursts/clusters…").

Notice how careful I am about proposing this: it is a mathematically very convenient model, and it's pretty accurate to model a lot of naturally occurring behaviours, but it can badly jump in your face if you forget it's likely a simplifying behaviour.

Take this example (which I made up right now, but I think it's realistic enough to make my point):

  Our packets are 0.1 s long. Interference emissions are much shorter, 0.001 s; they either appear within a packet (leading to a broken packet), or they don't. We model the occurrence of the interference to be a Poisson process.

This works well if the interferers are, for example, independent sensors that send a short packet when someone enters a room. In essence, your model represents the reality of when interference bursts happen accurately.

Even when the interferers are actually doing ALOHA for medium access, due to the packets being so short, the chance that you'll see them clustering is so low that you can still assume disjoint time periods (i.e. your packet transmissions) have independent interference arrival probabilities.

Now, if you happen to make your packets shorter, e.g. 0.008 s, you'll start "seeing" that ALOHA of the interferers mean that if there is indeed a correlation between interference occurrence during a packet and during a following packet.

If you made assumptions for your NOMA system based on the Poisson-ness of the process (like, you're doing CSMA/CD with deterministic backoff, you don't use an overly long spreader to distribute your info,…), you will run into a real congestion problem although your channel isn't actually heavily occupied.

So, treat your model assumptions with care. If in doubt, back them up with a quick simulation of what happens if they aren't fully met; if you can show that your overall system is robust about small model deviations, that greatly increases the credibility of whatever you're writing/designing.


Can I apply Fourier Transform to a non-time-indexed signal? Say I have a signal that is not x-indexed. That is, the x-axis of the signal is the distance traversed by car and the y-axis is the heading direction of the car at the corresponding distance.

Can I apply the Fourier Transform to this signal?
If so, what is the physical meaning of this transformation? I believe that the horizontal axis is no longer frequency any more. What is it in this case?

Yes you can. The unit of the "frequency" axis after the transform will be $m^{-1}$, and is known as spatial frequency.

For example, if there is a strong peak in the Fourier transform at $3.10^{-4} m^{-1}$, it means that your original curve exhibits a strong pattern that repeats at a scale of every $3.3km$, and from that you could infer that maybe the signal was recorded from a vehicle doing laps at the Monaco Grand Prix. The harmonics of this spatial frequency would contain a "signature" of the shape of the circuit.

A practical application of this is handwriting recognition - looking at shapes in the Fourier domain yields representations invariant to scaling, rotations, or more robust to deformations than the original data.


determining type of filter given its pole zero plot How can I classify a filter given its pole-zero map. For example I've got my zero's located at $\pm j$ and my poles located at $\pm\frac{1}{2}j$.

You'd have to figure out the frequency response of the filter. Here are two methods. I prefer Method 2 because it's quick and dirty, and you don't really care about the exact gain values in the frequency response, just the general shape to figure out the type of the filter.

Method 1: Brute Force/Computer Assisted

import scipy.signal as sp
import numpy as np
import matplotlib.pyplot as plt

b = [1,0,1]
a = [1,0,0.25]
w,h = sp.freqz(b,a)
plt.plot( w/np.pi, abs(h) )

Clearly there's a stop band around $\pi/2$, so this is a band-stop filter.

Method 2: Hand Computation

The frequency response can be obtained by moving along the unit circle in $z$-domain. The amplitude of the frequency response can be calculated by dividing the product of the distances from the zeros by the product of the distances from the poles.

In your example the pole-zero plot looks like this:

You can trace out the amplitude of the frequency response by moving the point $P$ along the unit circle from $0$ to $\pi$ and computing $\frac{\mbox{length(PO)}}{\mbox{length(PX)}}.$

Initially,  $\mathrm{length}(PO)>\mathrm{length}(PX)$, so the amplitude response is larger than $1$. Then gradually as we move closer to $\pi/2$, $\mathrm{length}(PO)$ keeps getting shorter until it becomes exactly zero at $\pi/2$.


Applying a high-pass filter on constant signal I recently implemented a Butterworth high-pass filter (2nd order). Everything seems to work fine, except that I have a question when applying such a filter to a constant signal. As there are no high-frequency components in a constant signal, I would expect the filter to yield a constant 0 signal.

In the plots below are my results. There seems to be ``ripple'' in the first couple of frames before the high-pass filtered signal (correctly) converges to 0.

Is this a logical result from applying a Butterworth high-pass filter on a constant signal, or might there be a bug in my code?

EDIT: I've created the same filter in Octave, resulting in the same output. This indicates that we're indeed looking at the step response, as Paul R. indicates.

Octave output

The input to your filter is essentially a step function (since it has value 0 prior to t = 0, and a positive value for t > 0), so you see the step response of the filter, hence the initial ringing. This is expected behaviour, and after a suitable amount of time the step response will have settled to zero.


Find continuous signal given a condition on its samples Let $x(t)$ be band-limited with $B = \omega_m$. Sampling gives us $$x(nT_s) = \begin{cases}
1,  & n = 0 \\
0, & n \not = 0
\end{cases}$$
And $\omega_s = 2\omega_m = \frac{2\pi}{T_s}$. Find signal $x(t)$.
My try: The first problem is about the definition for band-limited signal. It means $X(j\omega) = 0$ for $|\omega|\gt\omega_m$ or $X(j\omega) = 0$ for $|\omega|\ge\omega_m$? The sampling theorem requires that if $X(j\omega) = 0$ for $|\omega|\gt\omega_m$ then $\omega_s \gt 2\omega_m$ to avoid aliasing. So in this case is it possible to find interpolating functions other than $$x(t) = \begin{cases}
\frac{\sin(\omega_mt)}{\omega_mt},  & t \not= 0 \\
1, & t  = 0
\end{cases}$$which I've found using ideal low pass filter? I mean with the given information is $x(t)$ necessarily unique?

The condition
$$x(nT)=\delta[n]\tag{1}$$
is called the Nyquist criterion for zero intersymbol interference (ISI). It is important for the design of transmit pulses in digital communication systems.
Condition $(1)$ can be expressed in the frequency domain as
$$\frac{1}{T}\sum_{k=-\infty}^{\infty}X\left(\omega-\frac{2\pi k}{T}\right)=1\tag{2}$$
where $X(\omega)$ is the Fourier transform of $x(t)$. From $(2)$ we see that the shifted spectra need to add up to a constant. This is only possible if the maximum frequency of $x(t)$ satisfies $\omega_m\ge \pi/T$. For the minimum value $\omega_m=\pi/T$, the shifted spectra don't overlap and, consequently, $x(t)$ must be an ideal low pass signal with a flat spectrum. This is the only solution and it corresponds to the solution you came up with. If $\omega_m>\pi/T$, there are infinitely many solutions to $(2)$. One well-known example are the raised-cosine pulses.


How to Calculate the Period of the Discrete Time Sequence from Taking Its DFT This is my discrete time periodic signal, the time resolution of the samples is 0.2ms.

And it's periodic after each 11 samples. So the fundamental period of the signal is 2ms or the fundamental frequency of this sequence is 500Hz. The amplitude of each pulse increases from -5V with unit step to +5V.

But I want to find its frequency (which is 500 Hz) by calculating its FFT and plot of FFT spectrum. Please help me. Thank you.

Sample the signal into MATLAB.
Apply the fft function on it.

Plot the absolute value of the DFT signal (Its first half of samples).
Look for its peak, the index number of the peak multiplied by the sampling frequency divided by the number of samples will give the frequency.
Look at its inverse and your answer will pop.

If you upload a digital copy of the signal, I will show you in code.


How can I find the value of Inverse DFT I have dt signal $x[n]={[6.29, 8.11,-7.46,8.26,2.64,-8.04,-4.43,0.93,-9.29]}$

And I need to give the function value of:

1) sum of $x[n] =  \sum\limits_{k=0}^{N-1} X[k] \, e^{-j 3 \pi k/5}$ from k=0 to 9

How can I find that

I don't know if I understand exactly what you mean, because your question is very unclear (you may need to edit it), but I will try to answer what you may need in your task.

For DFT you will need this equation:

$ X[k] = \sum\limits_{n=0}^{N-1} x[n] \, e^{-j 2 \pi nk/N} $

where $k = 0, 1, . . . , N − 1$

But I guess you will be asking for IDFT more, so first you need to understand how IDFT works. I will try to explain this in my picture:

So the formula will look like this:

$ x[n] = \frac{1}{N} \sum\limits_{k=0}^{N-1} X[k] \, e^{j 2 \pi nk/N} $

where $n = 0, 1, . . . , N − 1$

Note that here, x[n] will be our complex time series, and X[k] is complex frequencies.

Now, we can get your sum, thanks to a simple code:

x = [6.29 8.11 -7.46 8.26 2.64 -8.04 -4.43 0.93 -9.15 9.29];
XF = fft(x);
k =  0:9;
YF = exp(-i*3*pi*k/5).*XF;
output = [sum(YF)];
disp(output)

And, as we can see, the sum equals 9.3.


How does parallel structure implementation of IIR speed up the process ? What happens in case of FIR? I was learning about digital filter structures and there I came to know that parallel implementation of IIR filters speeds up the process and I understand that.But what I am not convinced of is that parallel implementation of FIR filters doesn't speed up the process. I want to know the importance of parallel structures in FIR filters.
EDIT: I understand that if I have a digital IIR filter whose transfer function is $H(z)$ then its parallel realization would be
$$H(z) = \sum_{i}{H_{i}(z)}$$
where each $H_{i}(z)$ is no more complicated than bi-quadratic. So at the most we require two delay elements for each parallel part. So if order of $H(z)$ is 10 then ignoring time required for  multiplication, the minimum processing time is twice of the sampling period not ten times of that.
Now my question is what happens in case of digital FIR filters ?

As already pointed out in the comments, it's a misunderstanding that the parallel structure of IIR filters "speeds up the process". The reason for breaking up the transfer function into second-order (and first-order) blocks is to increase numerical stability, especially for fixed-point implementations. The low-order filter sections can be arranged in a cascade structure - i.e., effectively factoring the transfer function -, or in a parallel structure, which corresponds to a partial fraction expansion of the transfer function.

Note that if we don't take quantization effects into account, the direct implementation of the total transfer function, the cascade structure, and the parallel structure give exactly the same output signal. There is no difference in delay, as you seem to believe. Even though in the direct implementation you have a longer delay line, you don't need to wait till the input signal has traveled through the delay line, because the delay line is tapped. So the first non-zero input sample can directly produce a non-zero output sample.

Let's now come to your question about FIR filters. The transfer function of a causal $N^{th}$-order FIR filter with impulse response $h[n]$, $n=0,1,\ldots,N$ is given by

$$\begin{align}H(z)&=\frac{h[0]z^N+h[1]z^{N-1}+\ldots +h[N]}{z^N}\\&=h[0]+h[1]z^{-1}+\ldots +h[N]z^{-N}\end{align}\tag{1}$$

where the second line of Eq. $(1)$ is (trivially) a partial fraction expansion. That means that the "parallel structure" of an FIR filter simply equals a direct form implementation. In other words, there is no distinct parallel structure implementation of an FIR filter.


Origin of artifacts 12 Hz away from electrical line noise (50 Hz)? I made electroencephalographic measurements from people while they were sleeping. When plotting the frequencies, I see that I do not only have an expected 50 Hz (and 150 Hz and 250 Hz...) artifact (from nearby electrical devices) but also smaller artifacts to either side that are approximately 12 Hz apart (38 Hz and 62 Hz, 138 Hz and 162, ...) (examples from two participants):

Can someone tell me what might cause those artifacts?

Probably at exactly $\pm 12.5\text{ Hz}$, ie. some periodic phenomena that happen every four grid frequency cycles. Typical candidates for such things are rectifiers with little load.

At any rate, it's probably coming from your power grid, and you won't be able to get rid of them by anything short of notch filters.


Kalman Filter with Accelerometer with DC Offset Goal: For a particle moving uniaxially, to estimate position ($d$) and velocity ($v$) from noisy acceleration ($a$) and very noisy position (GPS) measurements using a Kalman filter.

Catch: The accelerometer has a DC offset, i.e. its zero is not zero but some number $\Omega$ such that the "true" acceleration $a^{true}$ is obtained from the measured acceleration coming out of the sensor $a^{obs}$ as:
$$
a^{true}=a^{obs}-\Omega
$$
My system states and transitions, note I've added the accelerometer bias or DC offset as a state are:

$$
\hat{x}_{k+1}=\left[\begin{array}{cccc}
d_{k+1} \\
v_{k+1} \\
a_{k+1}^{true} \\
\Omega
\end{array}\right]=\left[\begin{array}{cccc}
1 & dt & dt^2/2 & 0       \\
0 & 1 & dt & 0       \\
0 & 0 & 1 & 0       \\
0 & 0 & 0 & 1
\end{array}\right]\cdot\left[\begin{array}{cccc}
d_k \\
v_k \\
a_k^{true} \\
\Omega
\end{array}\right]
$$

And then the measurements are (notice I've used $a^{obs}$ instead of $a^{true}$):
$$
z_k=\left[\begin{array}{cccc}
1 & 0 & 0 & 0       \\
0 & 0 & 0 & 0       \\
0 & 0 & 1 & -1       \\
0 & 0 & 0 & 0
\end{array}\right]\cdot\left[\begin{array}{cccc}
d_k \\
0 \\
a_k^{obs} \\
0
\end{array}\right]
$$

Question 1: Does this look right and/or feasible and,

Question 2: How do I determine the process noise and measurement noise matrices $Q$ and $R$?

Thanks for your help.

Q1: Without actually implementing it and seeing that it works, your formulation looks pretty good... except that I do not understand your final equation:
$$
z_k=\left[\begin{array}{cccc}
1 & 0 & 0 & 0       \\
0 & 0 & 0 & 0       \\
0 & 0 & 1 & -1       \\
0 & 0 & 0 & 0
\end{array}\right]\cdot\left[\begin{array}{cccc}
d_k \\
0 \\
a_k^{obs} \\
0
\end{array}\right]
$$
That doesn't make sense because a) the output matrix is too big and b) $a_k^{obs}$ does not form part of the state, $a_k^{true}$ does. It should look more like:
$$
z_k=\left[\begin{array}{cccc}
1 & 0 & 0 & 0       \\
0 & 0 & 1 & +1       \\
\end{array}\right]\cdot\left[\begin{array}{cccc}
d_k \\
v_k \\
a_k^{true} \\
\Omega
\end{array}\right] = \left[\begin{array}{cccc}
d_k\\
a_k^{true} + \Omega
\end{array}\right]
$$
Notice the change in sign to get the expression for the acceleration measurement to come out correctly.

Q2: What noise is driving your system? Is it just the acceleration?

Your state update equation should look more like:

$$
\hat{x}_{k+1}=\left[\begin{array}{cccc}
1 & dt & dt^2/2 & 0       \\
0 & 1 & dt & 0       \\
0 & 0 & 1 & 0       \\
0 & 0 & 0 & 1
\end{array}\right]\cdot\left[\begin{array}{cccc}
d_k \\
v_k \\
a_k^{true} \\
\Omega
\end{array}\right]
+
Gw_k
$$
where $G$ determines which part of the state is driven by the noise process $w_k$.  The dimension of $w_k$ (is it a scalar or a vector?) will determine whether, if more than one part of the state is being driven, whether those drivers are different for each part.

If the driver is just acceleration, then choose:
$$
G = \left[\begin{array}{cccc}
0 \\
0 \\
1 \\
0
\end{array}\right]
$$
assuming that the offset in acceleration is not changing.

Then the question of what the covariance of $Q$ is boils down to: how quickly does your acceleration change?  That will really depend on your precise system. In this case, $w_k$ is just a scalar.

As to the covariance of $R$, that will depend on how different your measurements are form reality. I'd need more information to make a stab at that.


Down sample incoming stream of data to plot a dynamic graph I am working on a module where I need to plot ECG streaming data on a graph. This graph has configurable speed i.e. plot 1 heartbeat data sample within say 50mm,25mm, etc. For this, I am required to down sample based on the screen resolution and plot the same.

I am using Qt QML for plotting. What is the exact algorithm to achieve this? I don't have a clue on what exactly I have to google

Downsampling

Let's define our downsampling factor :
$$r = \frac{F_s'}{F_s}$$

Where $F_s$ is the original sampling frequency and $F_s'$ the reduced one.

If $F_s$ is a multiple of $F_s'$ (i.e. $1/r$ is an integer), then the dowsampling is straightforward : you just take one sample out of $1/r$. For example, if you want to half your sampling frequency, just discard every $2^{nd}$ sample and that's all (well not really, see next paragraph).

However, if that's not the case, you're going to have to interpolate between the samples. To make it simple, interpolating means trying to guess the value of a sample by knowing the value of the surrounding samples. You brain does interpolation all the time : when you see several points forming a line you instinctively draw the line in your head. Several types of interpolation exist, the more simple beeing the less precise. It usually depends on your application.

But there's more...

You've maybe heard of Shannon-Nyquist sampling theorem, which states that to sample a signal, you have to have a sampling frequency at least twice as big as your signal's highest harmonic. If you don't, it's going to cause something called aliasing. I won't explain this theorem here because there is plenty of information available everywhere.

But keeping this theorem in mind, you can probably understand that downsampling your signal might cause aliasing if you're not careful. The simple solution to that is to filter the signal prior to downsampling, at a cutoff frequency of $\frac{F_s'}{2}$. The steeper the filter's transition band, the less aliasing.

Back to your application

From what I understood, you have a signal that you want to map to a screen. For that, you need to downsample the signal to your screen's horizontal resolution. Let's call this resolution $H_x$ and your signal's sampling frequency $F_s$. For the sake of the example, let's fix $H_x=240\text{px}$ and $F_s=5\text{kHz}$.

You also have a time parameter (let's call it $t_r$) which sets the "zooming" of your scope, ie how long of a signal corresponds the length of the screen. Let's set it to $t_r=100\text{ms}$ (Disclaimer : I have absolutely no idea about the frequency of ECG so this is completely arbitrary).

That gives you :
$$
F_s' = \frac{H_x}{t_r} = \frac{240}{0.100}=2400\text{Hz}
$$
Note that that is not a multiple of $F_s$.

Applying what I said before, you have :
$$ r = \frac{F_s'}{F_s} = \frac{2400}{5000} = 0.48 $$

So first, you filter your signal with a brickwall at $f_c=\frac{F_s'}{2} = 1.2\text{kHz}$, then you downsample it by a factor of $0.48$ using interpolation.

Another thing to consider is the vertical resolution. For that it is a bit different as you are not downsampling the signal but quantizing it. However, the process is similar and I think you can extrapolate without our help.


complex numbers and fourier transform Is it possible to define a scaling property for fourier transform when the scale factor is complex?
Usually the scaling factor is real. What happen when a scaling factor is complex?

there are issues.  given this convention for the continuous Fourier transform (and inverse)

$$ \mathscr{F} \Big\{ x(t) \Big\} \triangleq X(f) \triangleq \int\limits_{-\infty}^{+\infty} x(t) e^{-j 2 \pi f t} \ \mathrm{d}t $$

$$ \mathscr{F}^{-1} \Big\{ X(f) \Big\} \triangleq x(t) = \int\limits_{-\infty}^{+\infty} X(f) e^{+j 2 \pi f t} \ \mathrm{d}f $$

it changes the path of integration from the real axis to something else.  this comes up when using this fact:

$$ \mathscr{F} \Big\{ e^{- \pi t^2} \Big\} = e^{- \pi f^2} $$

to get, along with using scaling, this result:

$$ \mathscr{F} \Big\{ e^{j \pi t^2} \Big\} = \sqrt{j} \, e^{-j \pi f^2} $$

which is a linearly-swept "chirp" and it's spectrum.


Is this jumbled image due to difference between Matlab and Python slicing? When I import an image from Matlab I get the following jumbled image

rather than the (224x224x3 RGB) image of a single elephant I expected

I write the image from Matlab as a binary file

function save_bin(filename, input_image)
  r = input_image(:,:,1); % dump image for input to python
  g = input_image(:,:,2);
  b = input_image(:,:,3);
  fid = fopen(filename,'wb');
  fwrite(fid, r(:), 'float32');
  fwrite(fid, g(:), 'float32');
  fwrite(fid, b(:), 'float32');
  fclose(fid);
endfunction

And read it back into Python as follows

m_resized = np.reshape(np.fromfile('m_resized.bin', dtype=np.float32),(224,224,3))
image = Image.fromarray(np.uint8(m_resized), 'RGB')
image.show()

I've investigated the slicing differences creating a 3x3x3 tensor in Matlab and importing it and there's a fundamental difference in the slicing syntax between Matlab Python but this alone doesn't explain  what's going wrong

Any ideas?

Note: The reason I'm using a binary format rather than a jpg image is that I want to compare a Matlab CNN model to a Numpy implementation at each stage and using an image is an intelligible means of comparison rather than comparing lots of numbers.

The answer is that the order in which pixels are accessed in Matlab is different from that in Python

The problem is resolved by specifying the order in the numpy reshape function to be the same as that of Matlab using the order='F' flag

Where ‘F’ means to read / write the elements using Fortran-like index order, with the first index changing fastest, and the last index changing slowest

def display_image(input_image) :
  image = Image.fromarray(np.uint8(input_image), 'RGB')
  image.show()

m_resized = np.reshape(np.fromfile('../resnet50_matlab/m_resized.bin', dtype=np.float32),(224,224,3), order='F')

display_image(m_resized)

This results in the correct image being displayed as shown in the image below

Rather than the mangled image displayed using

m_mangled = np.reshape(np.fromfile('../resnet50_matlab/m_resized.bin', dtype=np.float32),(224,224,3))

display_image(m_mangled)


What is the underlying concept behind Bandpass sampling? Can you suggest me some books/webpages on Bandpass sampling?

I undestand that if the signal is restricted between $f_L$ and $f_H$, then the minimum bandwidth required is $2(f_H - f_L)$. But say the centre frequency is in GHz and the bandwidth is in MHz, then by this theory I would have the sampling frequency to be in a few MHz. So in case of GHz centre frequency signal(say 1 GHz), I collect the samples every microsecond. But the signal is changing every nanosecond. I did not understand this concept.

Yes, the signal is changing "every nanosecond" for any non DC signal.  But if the bandwidth is really limited to 0.001 of the carrier frequency, then those changes are very predictable, as the changes result in a carrier waveform that is a nearly perfect sinusoid of a known frequency.  The signal can't differ too much from that highly predictable sinusoidal waveform unless the bandwidth is increased.  Assuming the bandwidth doesn't increase, any changes to the near perfect sinusoid are so small that one only needs to check the waveform every few hundred cycles of the sinusoid to see if it has drifted a bit from the expected point on the expected cycle.  Thus under-sampling at somewhere above 0.002 cf can tell you all you need to know in addition to knowing that the bandwidth was absolutely constrained to 0.001 cf or below.

But the sampling jitter requirement has to be on the order of a very tiny fraction of a nanosecond, or else one will completely miss the expected point on the expected cycle of the 1 Ghz sinusoid.


How does "image synthesis" work? I'm interested in how does "image synthesis" in audio synthesisers such as Metasynth work.

Image synthesis is when a picture is read into a program and "interpreted" into sound.

But how does it work?

There was something like that in Aphex Twin – Windowlicker (1999 single), visible in the spectrogram with a logarithmic frequency axis:

According to the album's Wikipedia page, the image-to-audio conversion was made with Metasynth. I don't know the specifics of Metasynth, but a similar effect can be achieved by calculating the inverse discrete Fourier transform (IDFT) of each pixel column and concatenating the resulting time-domain frames. To make it sound smoother, the time-domain frames can be cross-faded and the phases of the frequency bins can be randomized to disperse the times when the phases of the frequencies match. The image can be distorted so that the vertical coordinate corresponds to a logarithmic frequency scale.

The problem could be stated as: How can I create an audio signal that produces a spectrogram that looks like the input image. To make a really good conversion, similar methods could be used as in time-scale/pitch modification. Some of those try to maintain phase coherence of frequencies (particularly of those that fall between frequency bins), across time-domain frames. Some of those try to maintain the time coherence of transients – things that are more concentrated in time domain than in frequency domain.


Transfer function of resonant filter with 2 poles, peak at $f_0 = 500\text{ Hz}$, and $\Delta f = 32\text{ Hz}$ This a contest question. I'd like some help because I can't find any materials related to this topic.

https://www.qconcursos.com/questoes-de-concursos/questao/ecd6c966-51

My english translation:

Find the transfer function of a resonant filter with 2 poles, peak on $f_0 = 500\text{ Hz}$, $\Delta f = 32\text{ Hz}$, and sampling frequency $f_s = 10 \text{ kHz}$.

-> a) $H(z) = \frac{0.062}{1-1.8831z^{-1}+0.09801z^{-2}}$

b) $H(z) = \frac{0.081}{3-2.641z^{-3}+0.09801z^{-2}}$

c) $H(z) = \frac{0.082}{1-1.8831z^{-1}+0.09801z^{-2}}$

d) $H(z) = \frac{0.762}{1-2.8831z^{-1}+0.09801z^{-2}}$

e) $H(z) = \frac{0.262}{2-1.8831z^{-1}+0.09801z^{-2}}$

Any hints?

You mistyped a few crucial numbers (the coefficients of $z^{-2}$). With the correct coefficients as shown in the link in your question, you quickly see that choice $b)$ can be discarded because it has a third order polynomial in the denominator.

For a peaking resonant filter we need two complex conjugate poles, so we can discard choices $c)$ and $d)$ because they have real-valued poles (due to the negative coefficient of $z^{-2}$). So we're left with $a)$ and $e)$, for which we need to compute the angles of the poles, because this is (approximately) where the frequency response peak will occur. Having obtained those angles we just need to check which one is closer to the angle corresponding to a frequency of $f_0=500\textrm{ Hz}$ (with $f_s=10\textrm{ kHz}$).

So the angle we're looking for is

$$\theta=2\pi\frac{500}{10000}=\frac{\pi}{10}$$

and it turns out that the poles of transfer function $a)$ have just that angle.


Doubt in K-L Transform equations I have a doubt, it may be silly. Please explain me.

In K-L Transform, For forward transform we have the below equation

For inverse transform we have the below equation

In reverse transform why the 'A' transpose is not common to both y and m(mean of x) as in Forward transform equation?


Magnitude of the Gradient in Frequency Domain I'm learning some basics of image processing. Recently I've read about image filtering and two-dimensional Fourier transform, because I'm preparing for exam. And I have one question I don't know answer for: is it possible to do filtering in frequency domain (by using Fourier transform), which in spatial domain is implemented by computing magnitude of gradient of the image?

I know that gradient operator is defined:
$$
\nabla f \equiv grad(f) = \begin{bmatrix}g_{x} \\ g_{y}\end{bmatrix}
$$
and the magnitude of gradient:
$$
M(x,y) = \sqrt{g_{x}^{2} + g_{y}^{2}}
$$
In spatial domain we can use gradient to find edges. High frequencies are responsible for edges. High frequency means the rate of change is high, thus gradient has high value.

But can we do similar thing in frequency domain? Has the magnitude of gradient operator any equivalent in frequency domain?

Remember that as part of computing the gradients of an image, you are really just convolving your image, with (two) spatial filters. This then gives you two 'gradient' images, which you can then combine.

For example, to compute the x-direction-gradient-image $G_x$, you would convolve your image with a filter $f_x$ that might look like this $\begin{bmatrix} 1 & 0 & -1\end{bmatrix}$. Similarly, to attain the y-direction-gradient-image,  $G_y$, you might convolve your image with a filter $f_y$ that looks like this $\begin{bmatrix} 1 \\ 0 \\ -1\end{bmatrix}$.

Then, your final gradient image would be $G = \sqrt{G_x^2 + G_y^2}$.

So far so good. However, you don't need to do the convolution pixel by pixel. Instead, you could simply take the 2D DFT of $f_x$, multiply it pixel by pixel with the 2D DFT of the image, and then take the 2D IDFT back into the spatial domain. This would then also give you $G_x$. (Do the same to find $G_y$).

So the answer is at least a partial yes, in that you can compute the x and y gradient images via the FFT algorithm, (an implementation of the DFT), and then simply treating each of them as an element in a vector, and finding the L2 norm.


Audio interface or External sound card? I want to create a microphone array that can be connected to my laptop. Previous discussions in stack exchange and other forums suggested me to get a USB audio interface. I would want to do a lot of real time signal processing on the sound arriving at the microphones (To be precise, I want to calculate the time delay of arrival of sound at the different microphones). But from what I saw, the audio interfaces are a bit expensive. But external sound cards are not! Here:

http://www.ebay.com/bhp/external-sound-card-for-laptop

Can I use an external sound card to receive input from multiple microphone (create an array, that is).

I will require 4 mics, but right now I need to do it with only 2 mics to test the algorithms. Sampling rate would be around 40Khz.

Sound card == audio interface

The only difference is that the name "audio interface" is usually used for pro audio equipment. On the other hand, "sound card" usually refers to consumer audio.

Currently there are no cheap multi-input audio interfaces. You may want to check out second hand market, but even there USB and 1394 interfaces are expensive.

The only cheap solution that comes to my mind is to use multiple stereo-input interfaces and synchronize theirs clocks using SPDIF. Then you can link them on Linux (in asoundrc) or Mac OS. I don't know if it's possible on Windows.

However, I have never tried synchronizing clocks using SPDIF myself.

Without synchronizing clocks you need to use resampling which is computationally expensive and can degrade quality. alsa_in and alsa_out tools for Linux/JACK can do this. But they're usable only if you don't need accurate timing as they don't do latency compensation.


Understanding audio downsampling code (JavaScript) and adapting to a different ratio I inherited responsibility for some audio downsampling code that I don't fully understand:

  // downsampling variables
  var filter = [
      -0.037935, -0.00089024, 0.040173, 0.019989, 0.0047792, -0.058675, -0.056487,
      -0.0040653, 0.14527, 0.26927, 0.33913, 0.26927, 0.14527, -0.0040653, -0.056487,
      -0.058675, 0.0047792, 0.019989, 0.040173, -0.00089024, -0.037935
    ],
    samplingRateRatio = this.options.sourceSampleRate / TARGET_SAMPLE_RATE,
    nOutputSamples = Math.floor((buffer.length - filter.length) / (samplingRateRatio)) + 1,
    outputBuffer = new Float32Array(nOutputSamples);

  for (i = 0; i + filter.length - 1 < buffer.length; i++) {
    offset = Math.round(samplingRateRatio * i);
    var sample = 0;
    for (var j = 0; j < filter.length; ++j) {
      sample += buffer[offset + j] * filter[j];
    }
    outputBuffer[i] = sample;
  }

The audio is always in $32$-bit floats (from the Web Audio API).  The source sample rate is typically, but perhaps not always $48\textrm{ kHz}$. The output is always $16\textrm{ kHz}$ at the moment, although I'm interested in providing an $8\textrm{ kHz}$ option which is part of why I want to understand this better.

The things I do understand:

Each output sample is basically a weighted average of $21$ input samples
Correspondingly, most input samples are partially represented in $21$ output samples
Somehow it's filtering out frequencies above $8\textrm{ kHz}$, and that that is necessary to avoid aliasing
The weights (filter array) make an interesting graph:

Questions I have:

Is this a well-known algorithm that I just don't recognize because I'm a noob? (Perhaps an FIR Filter?)
Which part of the code is the low-pass filter? The fact that it uses $21$ samples when I'd expect $3$-$4$ to suffice? Or is it the "weights" rather than simply averaging the $21$ samples? All of it?
How would I go about choosing the filter values if I didn't already have this list?
Is it significant that the sum of filter array is almost-but-not-quite $1$? ($0.98198732$)
Besides the TARGET_SAMPLE_RATE, what would I need to change for $8\textrm{ kHz}$ output? I assume the low-pass filter would need to go down to a $4\textrm{ kHz}$ max, but I don't know how to do that.

Answers to one or all questions are welcome.

The vector of coefficients that you gave implement a finite impulse response (FIR) filter with a cutoff of 8 kHz (for an input rate of 48 kHz). As you surmised, these coefficients specify weights that are applied to a sliding window of input samples to generate each output sample; this is the basic operating principle of an FIR filter.

An illustrative way to look at the filter's behavior is via its frequency response. This indicates what frequencies are passed through the filter unchanged and which are attenuated. Generally, if you want to downsample a signal, you first need to bandlimit it so that any frequency content above the post-downsampling Nyquist rate are sufficiently attenuated, to avoid aliasing as much as needed for a particular application.

jojek provided the frequency response in a comment above, which I will shamelessly reproduce here:

Some observations:

This shows a lowpass characteristic, with a band of low frequencies passed through the filter with little to no attenuation, followed by a transition region where the frequency response falls off. It then enters a stopband where all frequencies are significantly more attenuated than in the passband. There are no hard-and-fast objective terms that can be used in all cases to identify exactly where these regions start and end, but you can get an idea by looking at the plot.
This is not a particularly good filter in terms of stopband attenuation (it only attenuates frequencies above the desired band by ~20 dB at minimum), but better specifications would require more coefficients (and therefore more computations per output sample).

Since the filter only attenuates the stopband content by ~20 dB, this data will alias into the signal that you are left with after downsampling. Any content that was in the stopband region in the original signal will be aliased down to some other frequency in the downsampled signal, with its power level decreased by ~20 dB (or more, depending on exactly where in the stopband it is; notice that there are some narrow regions, often called nulls or zeros of the filter, where the filter generates very high attenuation).

This is an equiripple filter, in that the ripples in the stopband have an equal peak height (at just under -20 dB). Based on some quick playing around with Matlab's Filter Design and Analysis tool, it looks like this was designed for a passband from 0-7 kHz and a stopband starting at 9 kHz, with equal weights to the passband and stopband approximation error (these are parameters for the Remez exchange algorithm used to design equiripple filters).

If you wanted to enhance or duplicate this design, you would need to use a digital filter design procedure. As far as implementation goes, these range from pretty simple (e.g. the window method) to quite complex (Remez exchange as I described before).

However, there are a lot of tools already out there to do this for you, both in paid tools like MATLAB and in free open-source tools like SciPy. Here's an example of how you could use SciPy to design a filter for the 4 kHz output rate that you mentioned (you'll also need matplotlib installed to plot the result):

from scipy import signal
import numpy as np

# Sample rate of the input signal (Hz).
fs = 48e3
# Desired number of coefficients in the lowpass filter.
N = 21
# Passband edge (Hz).
fpass = 3e3
# Stopband edge (Hz).
fstop = 5e3
# Generate the desired filter coefficients. The arguments can be interpreted as follows:
#   - the desired length of the filter
#   - list of band edges; this describes a passband from 0 Hz to `fpass` Hz
#     and a stopband from `fstop` Hz to `fs/2` Hz
#   - a response of 1 in the passband is desired (0 dB attenuation), and
#     a response of 0 in the stopband is desired
#   - since we passed frequencies in Hz to the filter design routine, it needs
#     to know the value of the sample rate
h = signal.remez(N, [0, fpass, fstop, fs/2], [1, 0], Hz = fs)

# Normalize the filter to 0 dB gain at zero frequency (which is just equal to
# the sum of the filter coefficients).
h = h / np.sum(h)

# Calculate the frequency response of the filter, then calculate its magnitude
# response, measured in dB.
freq, response = signal.freqz(h)
ampl = 20*np.log10(np.abs(response))

import matplotlib.pyplot as plt
fig = plt.figure()
ax1 = fig.add_subplot(111)
ax1.plot(freq/(2*np.pi)*fs, ampl, 'b-')
plt.title('Magnitude response of designed filter')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude (dB)')
plt.grid()
plt.show()

Note that digital filters don't have any inherent knowledge of their sample rate. Instead, their behavior is normalized to whatever sample rate you push through the filter. For example, if you took your original filter and passed data into it at 24 kHz instead of 48 kHz, it would have a cutoff commensurate with a 4 kHz output rate instead of the 8 kHz of the original.

With that said, I'm not saying that this is a great filter design. There is noticeable ripple in the passband and the stopband attenuation is relatively poor. There are a lot of parameters that you can adjust with this method, and there are lots of other filter design techniques as well; this should get you a start. As far as a text goes, I would recommend Understanding DSP by Lyons; it's a good introductory book.


How to make this cool Audio Spectrum Visualizer I want to make/create this spectrum visualization on a physical Led_Bar_Panel(s).
Here's what I think I should do:

Take sample's in real-time( Sampling freq = Fs KHz, only voice)
Do a FFT( As far I know we take (2 ^ n) samples at a time. Is there a mathematical relation between n and Fs ?)
Take the averaged amplitudes(eg. 1-10 Hz will constitute one band/bar of LED) amplitudes.
Pass amplitude information to an Arduino.
Arduino Controls the Bar lights as per incoming data from PC.

Is this approach right?

One other thing that I want to know is "Can a Arduino do this fft in real time @ 8Khz sampling"?

Kindly provide me information/suggestions.
Thanks in advance!

Is there a mathematical relation between n and Fs ?)

$n$ will determine your frequency resolution. For example, if $f_s = 16kHz$ and $2^n = 8192$, each bin in your FFT will have a width of about $16000 / 8192 \sim 2$ Hz. From your image, it looks like the smallest FFT bin width you want is 10 Hz wide. This means that if $f_s = 16kHz$, you can work with a size as low as $2^n = 2048$, but no smaller.

  One other thing that I want to know is "Can a Arduino do this fft in
  real time @ 8Khz sampling"?

There are many types of Arduino boards. This is absolutely doable with the recent boards that have 32-bit ARM processors. On the older 8-bit AVR-based Arduino boards, the data point I found is that a size 256 FFT runs in 7ms. So it would still be possible with some corner cutting, such as:

Reducing the sample-rate to 8kHz as you suggest, and get a few less bars in the upper range (or faking the data - there's not much going on in here which is not correlated with the few bars before it anyway...)
Reducing the FFT size and get a few less bars in the lower range
Reducing the frame rate at which you do the analysis, and if necessary use some interpolation between frames.

I'm pretty sure I've already seen a bunch of projects doing that. If you only need a few bars, another approach is to use a filter bank.


Why can't this system be adequately represented using a z-domain transfer function? According to this question and answer the following system cannot be adequately captured by a z-transform transfer function.
$$y[n] = y[n-1] + F_{\psi}(y[n-1)) + F_{\phi}(x[n-1])$$
where $F_{\alpha}(z)$ is a first order high pass filter of the form
$$F_{\alpha}(z) = \frac{\alpha (1 -z^{-1})}{1-\alpha z^{-1}} $$
The answer states that

The problem is that there's a pole-zero cancellation that snuck right by me and everyone else. It's evident in the left-hand side of (1), where the derivative of yk is the subject of the equation.

So the reason that you can't solve this problem as stated using the final value theorem is that you cannot adequately represent the system using a transfer function. There may be some way to save this within transfer function notation, but I just tried and failed at the first step, so I'm going to do it in state-space

What limitations of the z-transform (or others) require that this system be analyzed using alternative methods? What features of systems in general pose the same difficulty and why?

A transfer function describes an LTI system. As such, the given system can be described by a transfer function. However, if there are non-zero initial conditions, the system is no longer linear because there's a contribution in the output that does not depend on the input signal but only on the initial conditions. Consequently, the transfer function cannot directly be used to compute the system's response if there are non-zero initial conditions.
Nevertheless, the (unilateral) $\mathcal{Z}$-transform can still be used to compute the system's response, even with non-zero initial conditions by transforming the difference equation and using
$$\mathcal{Z}\big\{ y[n-k]\big\}=z^{-k}Y(z)+\sum_{m=0}^{k-1}z^{-m}y[m-k],\qquad k>0\tag{1}$$
EXAMPLE:
Let's use a simple example with a similar pole-zero cancellation as in the original problem to illustrate the point. Consider a system described by
$$y[n]-y[n-1]=\alpha \big(x[n]-x[n-1]\big)\tag{2}$$
The corresponding transfer function is
$$H(z)=\frac{Y(z)}{X(z)}=\frac{\alpha(1-z^{-1})}{1-z^{-1}}=\alpha\tag{3}$$
Clearly, $y[n]=\alpha x[n]$ is a solution of $(2)$. It is also the only solution if we require the system to be linear. However, it is not the only solution if we allow non-linear systems because there are infinitely many solutions of the form
$$y[n]=\alpha x[n]+c\tag{4}$$
with an arbitrary constant $c$. Note that these solutions cannot be inferred from the transfer function $(3)$.
Let's now use the $\mathcal{Z}$-transform to solve $(2)$ with initial conditions $y[-1]\neq 0$ and $x[-1]=0$. Transforming $(2)$ using $(1)$ gives
$$Y(z)(1-z^{-1})-y[-1]=\alpha X(z)(1-z^{-1})$$
which results in the following $\mathcal{Z}$-transform of the output:
$$Y(z)=\alpha X(z)+\frac{y[-1]}{1-z^{-1}}\tag{5}$$
In the time domain this becomes
$$y[n]=\alpha x[n]+y[-1]u[n]\tag{6}$$
where $u[n]$ is the unit step. Eq. $(6)$ is just a causal version of $(4)$.
This shows that the $\mathcal{Z}$-transform can be used to compute the system's response with non-zero initial conditions, even though the transfer function alone is inadequate for solving the problem.


Efficient FFT (or SFFT) for adjacent parts of signal with overlap? Imaging we have a signal $x$, which is segmented to 50% overlapped vectors $x_1,x_2,..,x_m$ , and we intent to compute FFT of each segment. Is there anyway that we can reduce computation of FFT of each block. I mean, if FFT of $x_1$ has taken and we know $x_2$ has a 50% overlap with $x_1$, so half of $x_2$, have already went through some FFT process.

Is there any method or modification that reuse this information to reduce computations of FFT of each segment based on information of the previous segment?

One method might be to split your length $N$ overlapping FFTs up, so that you compute the $\frac{N}{2}$-point FFT of each block of samples that come in. Then, in order to get the $N$-point FFT that you want, you can combine the last two $\frac{N}{2}$-point results with appropriate twiddle factors.

This would technically do what you want, but in most cases (especially on contemporary PC-like platforms with good FFT libraries), it's not going to be appreciably (if at all) faster than just doing the overlapping $N$-point FFTs.


generator matrix coefficient in convolutional code I cant determine a type of the following code:
G_1=1   // g_1=1
G_2=11  // g_2=x+1

Accorsing to description, it is convolutional code but I dont understand the type ( code rate)?

It is a systematic  rate-$\frac 12$ convolutional with constraint length $2$.


Bandwidth vs. sampling rate of DAC What is the difference between the sampling rate and the bandwidth of a DAC? I have found in one paper a DAC speced as 15GHz 3dB bandwidth, and 64GS/s sampling rate, and I haven't been able to relate them to each other in any way...

The bandwidth of a DAC is related to the analog filtering on the output. Most often (and in your example), this is specified as the 3 dB point of the analog filters that are part of the DAC. They are there to reject the Nyquist images of the signal that are part of the digital signal spectrum. Analog filters are not usually sharp enough to reject the first image you don't want while retaining the high frequencies of the baseband image that you might want to retain. For this reason, the 3 dB point starts a bit earlier to ensure that the required attenuation can be reached at the first folding band. Trying to output a signal with frequency components higher than this point will result in those components undergoing significant attenuation (resulting in a distorted output signal). Usually, one desires spectral flatness in the passband (low distortion of the time domain signal). In this case, it is recommended to oversample by at least 4 or 8 times (although the actual required amount is application specific). The bandwidth is important. In your example, if you are using a 64 GS/s DAC and want to use it to output signals higher than (or even close to) 15 GHz, then your signals will be output with significant attenuation.


How to select the order of LPC to model formants? I am having some difficulties understanding LPC and speech processing in general. The professor asked us to provide an heuristic way to select the order of LPC? For example if I wanted to model 3 formants what would be the order I should use? I understand the order of prediction is the number of coefficient used in the estimation of the signal $s(n)$:
$$\hat{s}(n)=\sum_{k=0}^p a_ks(n-k)$$
I am also open to any resources for a beginner in automatic speech recognition.

Usually you use 2 complex conjugate poles for each spectral peak. By other words, the order of the LPC model should be the double of the spectral peaks you want.
If you have 3 formants then you should use, at least, a 6th order LPC model.
You can try this yourself: record some vowels and then try (for each vowel) to use matlab or python to plot an LPC based spectrum. Change the order of the model to see if it can mimic the original spectrum.


signed 2's complement of a negative number for Integer = 3 and Fraction = 9 I am working on signed 2's complement for my fractions coefficients. Below code can be used for positive inputs. However for a negative 2's complement when I added LSB with 1, then an error pops. I think it is because of overflow. Now How can I convert negative coefficients using 2's complement for I and F of 3 and 9? and how to deal with overflow?
a = 0.7; % coefficient (actually for -0.7)
n = 3;         % number bits for integer part of your coefficient
m = 9;         % number bits for fraction part of your coefficient
% binary number
de2bi = fix(rem(a*pow2(-(n-1):m),2));
% signed 2's complement for negative numbers
negadd = [1];
negadd = [negadd, zeros(1, length(de2bi) - length(negadd))];
flipadd = fliplr(negadd);
magsig = fliplr(de2bi(bi2de(fliplr(de2bi))+bi2de(fliplr(flipadd))));
% back to binary
b2d = magsig*pow2(n-1:-1:-m).';

If your routine works correctly for converting positive numbers, you could just use it to convert a negative number by first adding the value of the sign bit (MSB) to the negative number (which must make it positive if it is not out of range), then convert the resulting positive number, and finally add the sign bit.
Example: convert the number $-0.7$ to Q3.9 format. The MSB has a value of $-2^2=-4$, so we first convert the number $-0.7+4=3.3$, which results in
011.010011010 ($= 3.30078125$)
Now you just add the sign bit and you obtain the representation of the original number:
111.010011010 ($=-0.69921875$)


Proof that first difference filter amplifies noise I'm a bit befuddled by noise's effect on derivative filters.

I've always 'known' that straightforward first difference derivative filters of discrete signals amplifies noise, but I'm struggling to prove that fact.

I thought it was caused by a greater gain for higher frequencies, but when I take the z-transform for a first difference filter I get
$$ H(Z) = 1 - Z^{-1} $$

and if I restrict the magnitude of $Z$ to 1 to get the DTFT and take $|H(Z)|$ the result is
$$ 2 |\sin(\frac{\omega}{2})|$$

which isn't increasing or constant in frequency like a high pass filter.

Could you explain why the noise increases if not for this effect?

Any filter whose frequency response magnitude is above unity at a given frequency $\omega$, amplifies signals at that frequency.

In discrete-time, the typical white noise is a wideband
signal whose power spectrum ranges from $-\pi$ to $\pi$.

The digital  first difference filter is an LTI system with the impulse response $h[n] = \delta[n] - \delta[n-1]$ and frequency response $H(\omega) = 1 - e^{-j\omega}$, whose magnitude is:

$$ |H(\omega)| = \sqrt{ 2 - 2 \cos(\omega) } = \sqrt{ 4 \sin^2(\omega/2) } = 2|\sin(\omega/2)|$$

as you can see as $\omega \to \pi$, then the magnitude goes to $2$ which is larger than unity and therefore amplifies any noise content found at those frequencies.

Note that, in any digital system, there will always be noise due to inherent quantization of numerical representation format, and therefore it will always be amplified using the given first difference filter.


Inverse Fourier transform of $\frac{\frac{2}{1+2i\pi f}-\frac{2}{4+2i\pi f}}{\frac{1}{1+2i\pi f}+\frac{1}{3+2i\pi f}}$ I wanted to calculated the inverse fourier transform of the transfer function :
\begin{align}
H(f) &= \frac{\frac{2}{1+2i\pi f}-\frac{2}{4+2i\pi f}}{\frac{1}{1+2i\pi f}+\frac{1}{3+2i\pi f}}\\
&\text{
I did some simplification :
}\\
&= \frac{3}{4}\left(\frac{1}{2+i\pi f}+\frac{1}{1+i\pi f}\right)
\end{align}
and after applying $\mathcal{F}^{-1}\{H(f)\}$ I obtained
$$ h(t)=\frac{3}{2}\left(e^{-4t}+e^{-2t}\right)u(t)$$ $$(this\ formula\ is\ corrected)$$
I used the fact that $$\mathcal{F}\left\{e^{-at} u(t)\right\}= \frac{1}{a+2i\pi f}$$
where $ u(t)$ is the unit step function.
Can I get some help to know if my process is correct?

This looks correct; obviously you need to simplify to obtain a transfer function with polynomials as numerator and denominator. Partial fraction expansion then leads to a form that can be directly transformed to the time domain. The only thing that might be wrong is the minus sign between the two exponentials.


Produce musical signal from instantaneous frequency I'm implementing the Yin Algorithm in Python so as to extract the pitch of my voice when whistling or humming. My goal is to produce the Gate (Envelop Control) and CV (Pitch Control) signals which can be inputted to an analog synth. I'd also like to write a software synth which takes theses signals as input. That way I can play the synth with my voice.
The top graph below is the output of the YIN algorithm when I whistled a few notes:

I then tidied this up with a low pass filter to produce the signal in the second graph. Finally and perhaps naively,  I produced the third signal as follows:
from scipy.io.wavfile import write
sr = 16000 # sample rate
f = [518.8087307738297, 518.9079938888592, 518.9177212823167, 519.8298830180304, 522.3027794382949, 523.8695096842832, 524.859458031283, 525.4703986196301, 525.4431662533539, 525.6859351990167, 525.6523745060124, 525.9697551477367, ...]
out = [math.cos(2*math.pi*f[i]*(i/sr)) for i in range(len(f))]

So in other words:
cos(2 pi f t)

where f(t) is the signal in the second graph. I expected this to produce a signal similar to that of my whistling recording put with constant amplitude. However, I get a really chirpy sounding signal that doesn't remind me at all of my original recording. If I use the first signal in my graph as f(t) I get noise when I apply it to cos(2 pi f t). What am I misunderstanding here? Here are the signals from samples 20000 to 30000. How should I reproduce a signal reminiscent of my original recording if I have the instantaneous frequency?

The f[i] * i together determine the phase and frequency at every i so sine wave is only rendered correctly if there is no speed change, so the algorithm must be changed to allow for continuous time i and any change in frequency must not cause a phase jump.
So try a phase accumulator approach, where phase for each i is incremented as required by the frequency.


What are the relative merits of various upsampling schemes? I recently encountered a DSP system which did some internal upsampling via zero padding.  Expecting zero-order-hold, I was surprised to find that a DC signal did not produce a DC output; many harmonics of the internal (lower) sampling frequency were also present in the output.

This leads to my question: What upsampling techniques are commonly used, and what are their relative merits?  Why would I choose zero padding, zero order hold, or first order hold, and what other techniques are available?

Some clarifications:

The system is real-time, so the upsampling scheme must be causal.
The upsampler is followed by a an anti-aliasing filter that can also be specified.

For the purposes of this answer I will use Matlab's terminology and define "upsampling" as the process of inserting $m-1$ zeros in between the input samples, and "interpolation" as the combined process of upsampling and filtering to remove the $m-1$ aliases ($m$ being the interpolation factor) that upsampling introduces.  For an explanation of how/why upsampling introduces aliases, please see this thread.

It is important to understand that any low-pass filter can be used to get rid of the aliases and thus complete the interpolation.  Some filters have advantages when used in interpolation, though.  I'll discuss the various flavors of interpolation filtering below.

FIR Filter

Interpolating FIR filters are efficient because they combine upsampling and alias filtering into one step.  This is most easily seen in an example.  Suppose we have a data sequence $x[n]$ and we want to interpolate it by a factor of two.  The first step is to upsample by a factor of two.  This changes the original data sequence from $x_0, x_1, ... x_N$ to $x_0, 0, x_1, 0, ... x_N$.

Now suppose we have a low-pass FIR filter, $h[n]$, that we will use to remove the alias.  When we convolve the upsampled data sequence with the filter, half the filter taps are stimulated by the non-zero samples, and half the taps are inactive because they correspond to the zero samples.  The half that is stimulated and the half that is inactive flips back and forth as the filter goes through the data.  These two sets of taps are sometimes referred to as the filter phases.

This same effect can be achieved implicitly by eliminating the upsampling, and filtering the original data sequence with an interpolating FIR filter.  The interpolating FIR filter produces $m$ outputs for every input sample.  For all $m$ outputs the filter will operate on the same $ceil(K/m)$ input samples (where K is the number of filter taps and "ceil" is the ceiling function).

An example will hopefully illustrate how this works.  Suppose that we have a six tap filter and we are interpolating by a factor of two.  The filter taps are [1 -2 4 4 -2 1].  If we literally interpolated and then filtered the samples and filter taps would line up (once there was full overlap) as follows:

$$
0 : 1\\
x_2 : -2 \\
0 : 4 \\
x_1 : 4 \\
0 : -2 \\
x_0 : 1 \\
$$
Next sample...

$$
x_3 : 1\\
0 : -2 \\
x_2 : 4 \\
0 : 4 \\
x_1 : -2 \\
0 : 1 \\
$$
Next sample...

$$
0 : 1\\
x_3 : -2 \\
0 : 4 \\
x_2 : 4 \\
0 : -2 \\
x_1 : 1 \\
$$
And so on.  The point of the interpolating filter is that it skips actually inserting the zeros and just alternates which set of taps it uses at the moment instead.  Thus, the preceding sequence would now look like the following:

$$
x_2 : -2 \\
x_1 : 4 \\
x_0 : 1 \\
$$

$$
x_3 : 1 \\
x_2 : 4 \\
x_1 : -2 \\
$$

$$
x_3 : -2 \\
x_2 : 4 \\
x_1 : 1 \\
$$

Zero Order Hold

A zero-order hold interpolator is one that simply repeats each sample $m-1$ times.  So a factor of two zero-order hold interpolator converts $x_0, x_1, ... x_N$ into $x_0, x_0, x_1, x_1, ... x_N, x_N$.  This method is attractive because it is exceedingly easy, both in terms of coding and computational load, to implement.

The problem with it is that its low-pass filtering is quite poor.  We can see that when we recognize that the zero-hold interpolator is a special case of FIR interpolation.  It corresponds to upsampling followed by an $m$-wide rectangle filter.  The Fourier transform of a rectangle filter is a sinc function, which is a rather shabby low-pass filter.  It's shabbiness can be fixed with a compensating FIR filter, but if you are going to do that you might as well just use a good low-pass filter to begin with.

First Order Hold

First order hold is a step up from the zero-hold interpolator in that it linearly interpolates the upsamples using the two nearest input samples.  So, a factor of two first-order hold interpolator would convert $x_0, x_1, ... x_N$ into $x_0, \frac{x_0+x_1}{2}, x_1, \frac{x_1+x_2}{2}, ... x_N$.

Like the zero-order hold interpolator, the first-order hold interpolator is a special case of FIR interpolation.  It corresponds to upsampling and filtering with a triangle filter.  For factor-of-two interpolation the filter is $[\frac{1}{2} 1 \frac{1}{2}]$, for factor-of-three interpolation the filter is $[\frac{1}{3} \frac{2}{3} 1 \frac{2}{3} \frac{1}{2}]$, and so on.

The triangle filter is two rectangle filters convolved together, which corresponds to sinc squared in the frequency domain.  This is a definite step up from the zero-order hold, but is still not great.

IIR Filter

I have never used an interpolating IIR filter so I won't say a lot about them.  I assume that the same arguments apply as in regular filtering- IIR filters are more efficient, can be unstable, don't have linear phase, etc.  I do not believe that they can combine the upsampling and filtering steps like a FIR filter can, but I could be wrong about that.

FFT Interpolation

I'll throw this one in even though it's not very common (of course, I don't think the zero-hold is common either).  This thread discusses FFT resampling, where resampling is both interpolation and decimation.

Higher Order Holds

Second-order hold interpolators are usually referred to as "quadratic interpolators".  They are non-linear and thus cannot be implemented as FIR filters, which are linear.  I do not understand the math behind them well, so I won't discuss their performance.  I will say, though, that I believe that they are somewhat common outside of signal processing.

Higher order (three or more) methods also exist.  These are referred to as "polynomial regressions".

EDIT:

Cascade Integrator Comb (CIC) Filters

I forgot to mention CIC Filters.  CIC filters are used for two reasons: they use only adders/subtracters (not as big a deal now that multiplies are fast and cheap), and they can do really large sample rate changes pretty efficiently.  Their down side is that they are essentially an efficient implementation of a cascaded rectangle filter, so they have all of the disadvantages of rectangle filters as discussed above.  CIC interpolators are pretty much always preceded by a compensating FIR filter that pre-distorts the signal to cancel out the distortion introduced by the CIC.  If the sample rate change is large enough the cost of the pre-distortion filter is worth it.


Why Convolve instead of Cross-Correlate for Detection of Known Signal? I am reading a journal on optimal signal detection by mobile sensor. It is mentioned:

  How should we plan to robots’ motion in order to maximize the detection probability? Assuming the source emits a known signal, the optimal detection algorithm is a matched filter (i.e., convolve the known waveform with the received signal and threshold).

Question: Why should we perform convolution of the known signal with the observed signal instead of cross-correlating the two, which shall detect the presence of similar patterns?

Convolving with Matched Filter is same as cross-correlation.

Suppose say your known signal is $x[n]$ for $0 \le n \le N-1$.

The matched filter is
$$
h[n] = x^*[N-n]
$$
where $*$ denotes conjugate (considering a generic complex signal). You can drop the conjugate for real signals.

The convolution with matched filtering operation is
$$
y[p] = \sum x[l]h[p-l]
$$
$$
h[p-l]=x^*[N-(p-l)]=x^*[l-p+N]
$$
Therefore
$$
y[p]=\sum x[l]x^*[l-p+N]=\sum x[l]x^*[l-\tau]
$$
where $\tau = p-N$ is the difference between index $p$ and $N$.

You can already see the matched filter output $y[p]=\sum x[l]x^*[l-\tau]$ is kind of cross correlation operation.

The Matched Filter output at $p=N$ is
$$
y[N]=\sum x[l]x^*[l-N+N]=\sum x[l]x^*[l-\tau]\Big|_{\tau=0}
$$
So as you can see here, Matched filtering output at $p=N$ which maximizes the SNR is the same as cross correlation output at zero lag, that is $\tau=0$.


Calculate Min/Max value of multiple sine waves added together So I feel like this should be possible, I'm just not sure how to do it. Let's say I'm generating a signal with a sample rate of 10K samples/second. Then let's say my signal is comprised of X number of frequencies( we'll use 3 for this example). For each sample I want to take the three frequencies and add them together, lets use 1000Hz, 1250Hz, and 1500Hz. Now that i've added them together, they might be too large (lets say the max value i can output is a 1) so then i need to normalize them. Is there a way to calculate the maximum or minimum value that these three frequencies will produce over a given period of time/number of samples so that i can get normalize my signal accordingly?

Thanks for any input on this!

The maximum value will occur when the frequencies all constructively interfere, which will certainly happen periodically if they are different frequencies.  Normalize by the sum of all the individual magnitudes.


How can I get the function of a curve from a dataset without using a curve fitting tool? previous relevant question here. I have a dataset of a signal captured from real hardware that produces a curve. Depending on the curve characteristics I can determine whether metal is present or not within a metal detector. The curve looks like this:

Each black line represent the average of a number of curves for each metal sample that goes through the system. The bottom curve being no metal present and the top, biggest curve being the largest metal piece present.
I want to remove the most significant portions of the curve and use definite integrals to help me decide whether metal is present. I can do this in matlab by putting the dataset of each average through the curve fitting tool, which generates the function (double exponential, 99.9% fit) of the curve and then from this I can put in my limits and integrate to get the area under the curve. I can then compare these areas to get good results.
My question: Is it possible to get the function of the curve without using a curve fitting technique?  Or is it possible to integrate the curve without curve fitting? I want to be able to analyse the signals with a (powerful) microcontroller in the end design and I do not want to write a curve fitting algorithm (time limited).

It's not necessary to first fit a curve to the data and then compute the integral. You can directly approximate the integral from the data using numerical integration methods. The most straightforward approximation would be a Riemann sum:
$$S=\sum_ny_n\Delta x\tag{1}$$
where $y_n$ are your data, and $\Delta x$ is the distance between the individual data points on the $x$-axis (assumed constant here, but that's not a necessity).
Of course there are more advanced methods of numerical integration, but $(1)$ is very simple and might be sufficient for your application.


pulse shaping filter and BER For example, in a PSK modulation. I saw that to plot the BER curves as a function of the noise, we never took into account the pulse shaping filter. More precisely, the signal is not filtered by the pulse shaping filter. In practice, it is necessary to filter the signal to have a limited bandwidth.
Why then is the pulse shaping filter not used to filter the sent signal in order to draw the BER curve?

Why then is the pulse shaping filter not used to filter the sent signal in order to draw the BER curve?

When drawing a BER curve, we're considering the bit or symbol energy.
That energy is the convolution of transmit pulse shape with receive filter evaluated at the symbol instance.
So, that filter is taken into account. We just usually define it such that the pulse energy is 1, and then you don't see it in your equations.


How can I use the fourier transform of a sequence of x-ray scan images to segment it? Sorry if this question seems to trivial, but I am a bit of a novice when it comes to signal and image processing and I need some guidance.

I have a 3D stack of about 256 grey scale images, each with pixel resolution 256x256, coming from an x-ray tomography scan of a rock sample.  I want to segment the data to isolate the void space from the solid space.  Typically I use ImageJ software to segment it using its 'auto' feature based upon a histogram of the grey-scale color values.  I am wondering if using fourier transform of the grey-scale data can somehow help me in the process.  What can the FT tell me about the image data that I can use for segmentation purposes?

Cited from wiki:

"a function (of time) is a representation of a signal with perfect time resolution, but no frequency information, while the Fourier transform has perfect frequency resolution, but no time information."

In your case, the Fourier transform has no spatial information. So for segmentation, you'd better try time-frequency analysis, such as wavelet transform


Obtaining beat histogram using Tzanetakis algorithm I am working to make a music player with automatic genre classification application in Java.
On Google searches I mostly found the works of George Tzanetakis who describe an algorithm from Music Genre Automatic Classification. I have big problems with obtaining the beat histogram.

http://dspace.library.uvic.ca:8080/bitstream/handle/1828/1344/tsap02gtzan.pdf?sequence=1  this is one of the documents that describes the algorithm and includes a flow diagram for beat analysis.

I have implemented all the functions from that flow diagram: the Full Wave Rectification, Low Pass Filtering, Downsampling, Normalization and Autocorrelation. I am pretty sure that these are correct because i found plenty of implementations and formulas in the internet.

What I find difficult is how to link these together to obtain the beat histogram.

What I do not figure is how to make the processing of those octaves and how to sum them. I understood that octaves are same thing as in music. I implemented the DWT transform that as it is described in this document and used the wikipedia page about DWT  implementation that luckily is written in Java, just like my application.

What I think my problem is  that I  did not understand properly how to get the octaves from the DWT transform. I obtain a waveform just like in : http://en.wikipedia.org/wiki/File:Haar_DWT_of_the_Sound_Waveform_%22I_Love_Wavelets%22.png

After that DWT transform I have split the transform just like in the wiki page with waveform for that audio sample: second half is first octave and just like in wiki example it looks like the original signal. First half  i split again and its second half is second octave and so on.   Am I doing wrong at this stage? For about two weeks I was sure I did it correctly and now that I get a wrong signal after autocorrelation, with no peaks as it should be.

How I did to sum the octaves after processing each of them as in the diagram: I first took the largest one,  and for the second one which is half size of first octave I added to the large one every second index in the array. From the third I added the values to every 4th index in the array and so on.  I got about 6-7 octaves in total. And this is where I think I may be wrong again. Because those octaves are arrays with different sizes(each is half size of the other one) and had to split the smaller ones over the larger octave. Have I done wrong? I  very much believe that the fact that my octaves arrays have different sizes could be where I am wrong.

Please tell me where I may be wrong in what I described about how I tried to implement that algorithm. I very much think that I did not understood what octaves  and I should not have arrays with different sizes after the DWT transform.

After the autocorrelation I should have a signal with peaks while I have something linear so I cannot make Peak Detection(last stage of the algorithm) with this.

One more thing I would like to ask is how would I get that bpm from the signal after the autocorrelation? I have an array with size of 2048  because I started with 131072 samples and did a downsampling to have 65536 values as in the documentation(my wav files have 44100 frequency while in document they had 22050). After dwt and getting first octave I had 32768 samples for first octave and 2048 after the downsampling by 16. In the beat histograms shown in this document they have bpm values from 60 to 200. How am I obtaining those values? I understand that bpm is almost same thing as frequency but I do not know how to get frequency from this transform.

Edit

After summing the octaves, using the sample from the answer, I got something like this:

They are a bit different because I used DWT as filterbank. After I apply the usual autocorrelation I get something like this:

But the problem is that Tzanetakis used the enhanced autocorrelation which is a bit different and is based on Tolonen and Karjalainen   algorithm.

The firsts steps by Tzanetakis can be described as an Scheirer algorithm(Tempo and beat analysis of acoustic musical signals by Eric Sheirer), you can use DTW, filterbanks or FFT to split your signal in N subbands!

I'll try to help you with what I did in the past, my steps are described below.

I can show a simple example of 6(six) (200, 400, 800, 1600, 3200, 6400) decomposed subbands aka "Scheirer 1998″, I used this sample to test:

After DTW or filterbanks you need extract the envelope of each subband:

Apply Full Wave Rectification
Apply Low Pass Filter
Apply Downsample
Apply Mean Removal

My envelopes from the six subbands above are:

OK now do you have six downsampled envelopes (Downsample can help you gain speed to compute the autocorrelation), you need sum all envelopes like:

Almost there, you need apply autocorrelation and to find bpm values between 60 and 200, you need compute the START and END position to be computeded from your summed envelopes:

start = (60 * downsampled_value) / (200)

end = (60 * downsampled_value) / (60)

Now you know the start and end to do the autocorrelation :-)

The next step is find the peak from your autocorrelation and convert to BPM:

BPM=60*downsampled_value/PEAK

And finally ride a histogram with the results of BPM's from each audio piece analyzed !

PS: Tzanetakis and Sheirer algorithm is not perfect its works well in strong rhythmic songs.


Why transforming the data to a high-dimensional feature space in which classes are linearly separable leads to overfitting? I read in my book (statistical pattern classification by Webb and Wiley) in the section about SVMs and linearly non-separable data:

  In many real-world practical problems there will be no linear boundary separating the
  classes and the problem of searching for an optimal separating hyperplane is meaningless.
  Even if we were to use sophisticated feature vectors, $\Phi(x)$, to transform the data to a
  high-dimensional feature space in which classes are linearly separable, this would lead
  to an over-fitting of the data and hence poor generalization ability.

Why transforming the data to a high-dimensional feature space  in which classes are linearly separable leads to overfitting and poor generalization ability?

@ffriend has a good post about it, but generally speaking, if you transform to a high dimensional feature space and train from there, the learning algorithm is 'forced' to take into account the higher-space features, even though they might have nothing to do with the original data, and offer no predictive qualities.

This means that you are not going to be properly generalizing a learning rule when training.

Take an intuitive example: Suppose you wanted to predict weight from height. You have all this data, corresponding to people's weights and heights. Let us say that very generally, they follow a linear relationship. That is, you can describe weight (W) and height (H) as:

$$
W = mH - b
$$

, where $m$ is the slope of your linear equation, and $b$ is the y-intercept, or in this case, the W-intercept.

Let us say that you are a seasoned biologist, and that you know that the relationship is linear. Your data looks like a scatter plot trending upwards. If you keep the data in the 2-dimensional space, you will fit a line through it. It might not hit all the points, but thats ok - you know that the relationship is linear, and you want a good approximation anyway.

Now lets say that you took this 2-dimensional data and transformed it into higher dimensional space. So instead of only $H$, you also add 5 more dimensions, $H^2$, $H^3$, $H^4$, $H^5$, and $\sqrt{H^2 + H^7}$.

Now you go and find co-efficients of the polynomial to fit this data. That is, you want to find co-efficients $c_i$ for the this polynomial that 'best fits' the data:

$$
W = c_1H + c_2H^2 + c_3H^3 + c_4H^4 + c_5H^5 + c_6\sqrt{H^2+H^7}
$$

If you do that, what kind of line would you get? You would get one that looked a lot like @ffriend 's far right plot. You have overfit the data, because you 'forced' your learning algorithm to take into account higher order polynomials that have nothing to do with anything. Biologically speaking, weight just depends on height linearly. It doesnt depend on $\sqrt{H^2 + H^7}$ or any higher order nonsense.

This is why if you transform the data to higher order dimensions blindly, you run a very bad risk of overfitting, and not generalizing.


How do I find the Energy Density Function of $g(t)$ if i am not given an input or impulse response? $$g(t)=\frac{12a}{t^2+a^2}$$
I need to find the Energy Density Function of the signal, but everywhere I look has an input and an impulse response. Does anyone know how to solve this. Would I just take the absolute value of it, square it and then integrate it?

If you're referring to energy spectral density, then you need to compute the squared magnitude of the Fourier transform of $g(t)$:
$$\big|G(\omega)\big|^2=\left|\int_{-\infty}^{\infty}g(t)e^{-j\omega t}dt\right|^2\tag{1}$$


Estimate changing time lag between signals Suppose I have two timeseries $x_1(t)$ and $x_2(t)$, shown in the image that I drew below. They look almost identical in general form, except some features are shifted slightly in time from one series to the other by an about $\delta t$. Sometimes they are measured first in $x_1$, and sometimes in $x_2$, so sometimes $\delta t$ is positive and sometimes negative. The change in lag between features is random, but constrained to be of a small range so that shifts never overlap and cancel.

I want to know, what is the mean time lag between the two signals?

Even better, would be to know the distribution of time lags so I could estimate the mean, variance and range.

Does anyone have ideas or suggestions for how to do this?

My own thoughts so far are: I can imagine calculating the cross correlation between signals $x_1$ and $x_2$. However, I would expect the result to be strong correlation at zero time lag, because the time lags can be both positive and negative, so that the average lag is zero. However, the width/stretching of the zero peak, must be related to the average time lag, however, I can't figure out how exactly at the moment.

Another thought I have had is to create a moving correlation estimator. In a similar way to how wavelet analysis can be used to estimate how the power of certain frequency components changes as a function of time. I could develop a function that measures the lag correlation as a function of time by breaking one of the signals into smaller chunks and convolving it across the second signal while stepping along.

All suggestions appreciated.

If the signals are as you've drawn them (flat, with abrupt changes), then perhaps this approach might work:

Take the absolute value of the differences between successive time samples.
The non-zero entries are the transition times.
Find these times for both signals.
Subtract the two.

I've made an attempt to do this in the R code below.

Example output:

 > print(differences_in_delay)
  [1]  2 -5  2 -4 -1  5 -4  3 -3  0
 > print(t[which(abs(diff(x2))>0)] - t[which(abs(diff(x1))>0)])
  [1]  2 -5  2 -4 -1  5 -4  3 -3  0

It breaks down when the transitions are within a sample or two of each other, but otherwise seems to work OK.

Will also probably break down if there's noise on your signals.

R Code Below

T <- 1024
t <- seq(1,T,1)
Nchanges <- 10
times_for_changes_for_x1 <- sort(runif(Nchanges,T/10,9*T/10))
differences_in_delay <- round(runif(Nchanges,-5,5))
times_for_changes_for_x2 <- times_for_changes_for_x1 + differences_in_delay
changes_for_x1 <- runif(Nchanges,-5,5)

x1 <- 0*t
x2 <- 0*t
change_no_x1 <- 1
change_no_x2 <- 1

for (tm in t)
{
  if ( (change_no_x1 <= Nchanges) && (tm > times_for_changes_for_x1[change_no_x1]))
  {
    x1[tm] <- x1[tm] + changes_for_x1[change_no_x1]
    change_no_x1 <- change_no_x1 + 1
  }
  else
  {
    if (tm > 1)
    {
      x1[tm] = x1[tm-1]
    }
  }
  if ( (change_no_x2 <= Nchanges) && (tm > times_for_changes_for_x2[change_no_x2]))
  {
    x2[tm] <- x2[tm] + changes_for_x1[change_no_x2]
    change_no_x2 <- change_no_x2 + 1
  }
  else
  {
    if (tm > 1)
    {
      x2[tm] = x2[tm-1]
    }
  }
}

plot(x1,type='l')
lines(x2,col="red")

plot(abs(diff(x1)>0)*1, type='l')
lines(abs(diff(x2)>0)*1, col="red")

print(differences_in_delay)
print(t[which(abs(diff(x2))>0)] - t[which(abs(diff(x1))>0)])


Effect of rounding off ADC value I have a 10-bit ADC with 8.5 bit ENOB. I have 2.5 quantization levels of thermal noise in front of the ADC. Is there a reason to process more than 8 bits from this ADC?

How wide is the word width of the processor receiving these 10-bit values from the A/D converter?  Wider than 8 bits?

What is there to gain by adding another layer of quantization error?

What you have here is a highly dithered signal that is dithered naturally.  Sometimes in audio we actually go out of our way to dither a signal so that further processing and quantization is better.

If your processor is 8 bits wide and you absolutely do not want to do multi-byte arithmetic (with an Add-with-carry instruction) or do not have the computational bandwidth to do higher precision arithmetic, then perhaps it makes sense to quantize your noisy 10-bit value to an 8-bit value.  Otherwise it does not.


Why did my sine wave turn into a square wave when written to a WAV file in Octave? I am trying to use Octave to generate a pure sine wave. The code for the same is as follows:

x = 10.*sin(2*pi*(300/16000)*(0:1:400));
The sampling rate is 16000Hz, the sine wave is at 300 Hz. I write the above wave to a file using wavwrite like so:
wavwrite(x, 16000, 16, "temp.wav")

When I try to read it back into a variable, like so:
y = wavread('temp.wav');, I get square waves upon plotting y.
I have checked the sine wave and the period indicates a frequency of 300Hz.

How can a pure sine wave become a square wave on simply writing and reading? Or am I going wrong somewhere?

x is being clipped and that is why y looks like a square wave. When writing x to disk using wavwrite, the samples of x are stored in 16 bits Q15 fixed-point format. That means your data must be in range -1 to +1 (in principle +1 minus one lsb). Therefore, x must be normalized to be in this range before calling wavwrite in order to avoid clipping.


Predistortion of signal constellation? I am trying to learn aspects of space/satellite communications. So Amplitude-Phase Shift Keying is quite popular in this domain.

APSK over Non Linear channels has worse performance than the earlier used PSK. Which is important because to non linear effects caused by the power amplifier on board. To counteract the effects of non linear distortion, one of the techniques used is Pre Compensation at the modulator.

This involves constellation predistortion. What does this mean?

Does this refer to error control encoding schemes like convolutional coding, tcm, turbo coding etc?

In very general terms, pre-distortion means altering a signal's shape, so that after it goes through the channel, the output looks like the original signal. So you basically alter the signal by the inverse of the distortion.

This is a separate issue from coding.


Measuring walking speed with iPhone accelerometer As part of a research project, I have been tasked with determining gait speed (walking speed) for a subject using a handheld iPhone. The main constraint is that walking distances are very short (10-30 meters), and accuracy is critical (data is used for research.) Using a pedometer and measuring step distance is not an option (too time consuming.) The only plausible way I foresee this could be work is using accelerometer data. However, I'm not sure whether this will succeed. I was wondering if anyone has experience with these sensors and could answer the following two questions:

Will an iPhone accelerometer plausibly be able to accurately measure walking distance for such short distances? Are there any figures on accelerometer accuracy for this kind of sensor? (I couldn't find any.)
Are there any established protocols for converting accelerometer data into distance measurements (calibration, noise correction, sampling rates, etc.) that have been shown to have good accuracy for such short distances? I have been able to find this, http://www.freescale.com/files/sensors/doc/app_note/AN3397.pdf but I'm not sure if this will apply well to the iPhone sensors.

I don't have any idea about Accelerometer in iPhone, but converting Accelerometer data into distance is fairly difficult because the offset in double integration causes drift as already told. Check this interesting presentation sensor fusion for further details.
As far as I  know this can be avoided by using Digital Accelerometers( with high sampling rates), double integration has to be performed using electronic integrators with appropriate filters at every stage to remove the offset as well as the low frequency noise.
Also in that presentation he mentions the possibility of using more sophisticated algorithms like use of Kalman Filter for object
tracking or Pedometer Algorithm for pedometer design.

Hope this gives little more insight into your problem:)

PS: This does not accurately address your question still I thought I will share my knowledge, I want to comment but I don't have enough reputation !


Bright lights on colors If you shine a bright enough white light on a red object, will the object turn white or stay red? Please assume the object is not capable of specular reflection.

When an object is Red, we mean that it only reflects the portion of the electromagnetic spectrum representing this color, about 650nm.

For a red object, we generally consider that only the color red is reflected, like the XKCD-like diagram below.

Source:

Color and Color Vision
Quora

EDIT:

This is a wavelength-intensity graph for various colored objects (anodised aluminum to be exact). The intensity is maximum for the wavelength associated with the color we see, but it still can reflect other colors, although at a much lower intensity as well.

Source:

Chen, Yuting, et al. "On the generation of interferometric colors in high purity and technical grade aluminum: an alternative green process for metal finishing industry." Electrochimica Acta 174 (2015): 672-681.


Transforming RGB images including NIR to LAB I'm relatively new to image processing, so I hope I don't ask trivial questions.

I have some images that I want to use in a machine learning context. The images have four color channels: RGB and NIR (near infrared). After some research, it seems like the LAB color space is better suited for machine learning in many cases than the RGB one.

First question would be: is that correct, or did I just fall for single people claiming something?

The second, more important question: Is there a way to properly include NIR in the LAB color space? Or would it be better to keep it as separate information?

And lastly, another question that came up: Can I even transform RGB to LAB properly? What information would I need (wavelengths of red, green and blue; some information on the whitepoint, ...)? Is there a different color space better suited?

If relevant: I am focusing on plant images and tasks appearing there, e.g. plant/background segmentation, fruit detection, leave counting,...

Let's start with the source of your RGB + near IR images. To break the process down:

Some camera must have taken those images.
The camera sensor is sensitive to light within two given wavelengths. The technical term for this is spectral response.
Alternating elements (Bayer pattern) in a sensor are etched with filters. Each filter type passes certain wavelengths while supressess the others.
The alternating sensor elements accumulate the filtered light from different wavelengths. This creates our very first raw signal. This raw space is then mapped onto a color space, most commonly the sRGB space.
Filter properties are precisely measured by the manufacturers. They use those properties (and assume a D65 white point, unless manually chosen to something else) to convert the captured raw space into sRGB which gets passed down to you (usually with a little compression and anti-aliasing filtering).

Answering your final question first, with the sRGB assumption in place, you are now armed with one very important source of information - your white point is D65 unless explicitly stated elsewhere. You may now safely proceed to convert your sRGB image into LAB. MATLAB has a handy function for the same.

To answer your questions:

  First question would be: is that correct, or did I just fall for
  single people claiming something?

A: People say a lot of things - largely from their personal experience. You must create your own. I personally exect LAB to help. You can read more about that on my Quora answer here.

  Is there a way to properly include NIR in the LAB color space? Or
  would it be better to keep it as separate information?

All traditional color spaces are designed to fall within the locus of human perception. Since near IR is out of bounds for our perception, it's usually modelled as a standalone entity (seperate channel). However, many imaging devices have a near IR spectral response below 750 nm. If that's the case, I would expect statistical correlation between your G or R channels with near IR. A simple PCA would answer that question.