Why this DTFT plot from CTFT? I have a sample test with an answer but don't understand how they got to the answer:

  $x(t)$ has info only between $2 < |\omega| < 4$
  $X^F(\omega) = 0$
  for other frequencies.

  All noise was cleaned from the signal and its CTFT is defined:
  $X^F(\omega) = \begin{cases} |\omega|-2  & \text{for 2 < |$\omega$| < 4}\\ 0 & \text{otherwise} \end{cases}$

  The clean signal was sampled
  at $\frac 1 T = \frac {12}\pi$
  Draw the DTFT of $X^f(\theta) \text{ for} -3\pi < \theta < 3\pi$

  Use a clear coordinate system.

The answer has x-axis span of -9 to 9 and y-axis of 0 to 4.

How was this answer calculated?

I thought I would get a simple repeating copy of the original CTFT, but the original CTFT with $2 < |w| < 4$ is a line beginning at w=2 (on the x-axis) with the value zero, going up until the point  $w=4$ where it has the value of 2...

Please help!!  Thanks

When a continuous time signal is sampled in discrete intervals, the effect in the frequency domain is that the CTFT becomes periodic in frequency - that is, you end up with "copies" of the original CTFT, with each copy centered at frequencies equal to k*1/T (k = any integer).

So, when you sample at 12/pi you will get a copy of the original CTFT centered at w=0, w=12/pi, w=-12/pi, w=24/pi, w=-24/pi, and so on.

Since you have frequency content between 2<|w|<4 in the w=0 case, the w=k*12/pi cases will have frequency content between 2<|w-k*12/pi|<4.

Be careful - if frequency content from multiple copies have frequency content in the same places on the frequency axis, the effect is that they add together.  This is called "aliasing" and is usually an unwanted effect of sampling.

Here is a nice graph that shows the effect of sampling a continuous time signal on the FT.

Thanks @rothloup !
I'm adding a layman's explanation to your answer and explaining the result:

There are four types of Fourier transforms: (that's easy to remember!)

All Fourier transforms work with infinite values so either we add zeroes around the given signal or we repeat the signal to infinity, with the length of the given signal as its periodicity.

The four types of transforms are:

CTFT or simply the FT (Fourier Transform): Continuous, Aperiodic
Used for a continuous signal with "one time" information in a small region of time. Zeroes added before and after to infinity, for using this transformation.
FS or Fourier Series: Continuous, Periodic
Used for a continuous periodic signal. We repeat the original (limited time region) signal before and after to infinity, for using this transformation.
DTFT: Discrete, Aperiodic
Used for a discrete signal with "one time" information in a small region of time. Zeroes added before and after to infinity, for using this transformation.
DFT: Discrete, Periodic
Used for a periodic discrete signal. We repeat the original (limited time region) signal before and after to infinity, for using this transformation.

The input of the DTFT - is discrete (sampled) input from an aperiodic signal (actually: information in a specified finite region of time, the rest set to zero...).

Because of the infinite zeroes being added, the period becomes extremely large, actually infinitely large, so there is only one period - or in other words the input is actually aperiodic!

Because of the zeroes added, the spacing in the frequency domain of the transform's output become extremely (actually: infinitely) close to each other.

So the output of the DTFT is Continuous and Periodic and also Symetric around the Nyquist number (in our case zero).

(That is why in the computer, where you deal only with sampled discrete data, we use the DFT and not the DTFT)

So if we have the CTFT $X[(\omega]$, we get the DTFT $X(\Omega)$ by setting:
$$\Omega = \omega T$$

In our case: since the samples were $12/\pi$ (in words: 12 times each cycle) and since the DTFT is symetric the period is half of the cycle (the other half is simply its mirror image) so the period is 6 and we multiply the $\Omega$ by that. While the areas with valid info are reduced: Formerly from 2-4, a quarter of the -4 to 4 period span, And become a sixth of the now -3 to 3 span:

$$ X(\Omega) = \begin{cases} 6|\Omega|-2  & \text{for 1 < |$\Omega$| < 2}\\ 0 & \text{otherwise} \end{cases}$$

a line from 4 to 0 between $\Omega$ 1 and 2
a line from 0 to 4 between $\omega$ -2 to -1. (actually down from 4 to 0 going backwards, from $\Omega$ 2 to 1)
and zero: between 0 and 1, between 2 and 3, and between 0 to -1 and -2 to -3.

As seen on the plot between -3 and +3

The pattern is repeated, starting with -3, resulting in:

$$ X(\Omega) = \begin{cases} \frac{6|\Omega|}{k}-2  & \text{for: $6k+1$ < |$\Omega$| < $6k+2$ ($k$ being any positive integer)} \\ 0 & \text{otherwise} \end{cases}$$

As plotted in the answer


How should calibration points be distributed? I do camera calibration for a single camera with a chessboard pattern in opencv for the first time.
The focus is fixed.
I'm not sure how to do it, there seems to be many questions:

How many images should be used for calibration?
What pattern should be used? number of points, size of squares
What transformations should be applied to the calibration target? rotated, tilted, shifted...

I came to the conclusion that none of the questions above can be answered on their own, because they all depend on each other. Say for example taking a lot of images that are very similar to each other are intuitively not as good as having less images that are better distributed in the frame.

I took a step back and thought about what's going on. (I tried)
Basically, there's a mathematical model of the camera and its aberrations (focal length, distortion, etc...) and the calibration is the process of collecting data to fit the model (its values) to the hardware at hand. This appears to be mostly statistics. In order to fit the model to the entire image well, sample data should be collected from the entire image.

I conclude that all questions above boil down to how should the detectable points of the pattern be distributed in the 3D space.

I'd like to know if this thought makes any sense.

I'm tired of reading handwavy suggestions like "the more images, the better", "at least arbitrary number of images should be used", etc. that do not provide a quantifiable measure for the quality of a calibration. Sure, there's a reprojection error, that tells how good the calibration fits to the provided calibration data, but not how good that data is suited for calibration in the first place.

I came to the conclusion that it makes more sense to evaluate the quality of a calibration by the distribution of the sample data used to create it. Makes sense?

Questions 1 and 2 above are thus combined to : What should be the overall density (or density distribution) of detected calibration points among all calibration images?

Say I take all detected points from all calibration images and place them in one image, I could see how many calibration points there are per pixel area $\text{px}^{2}$.If I have a lot of images with the pattern in the center of the image, I get a higher density there.

What would be a good minimum (baseline) value for the density of calibration points per pixel area?
$0.1 \ \text{px}^{-2}$ ? $0.5 \ \text{px}^{-2}$ ?

Should the distribution of calibration points per pixel area be constant over the entire image? Or do I need a higher density (read: more sample values) at the edges for example?

Question 3 is asking for the distribution in the 3rd dimension. The calibration target with the pattern can be moved back and forth (along the optical axis) or tilted. Both transformations lead to calibration points in front or behind the plane that's in focus1. Depending on the depth of field, corner points can be detected up to some distance from the focal plane. This adds a 3rd dimension to the space in which calibration points can be detected.
How should calibration points be distributed with respect to distance from the plane of focus1? Should most points be around the plane of focus1 given that detection of points further away from it is not as good?

tl, dr;

In front of a camera, there's a 3D space enclosed by a frustum of a pyramid/cone in which calibration points can be detected. How should calibration points be distributed in that space in order to get a good calibration?

1 With that I mean the plane of points that will be "in focus" in the image.

For that I have some guidelines, which I'd like to share. Even though some of these might not be valid for certain calibration software, for typical cases, they hold:

1) The entire pattern should be visible in all sequences.
2) At least 1/4th of the image should be occupied by the pattern.
3) Preferably, the lighting should be uniform everywhere on the pattern.
(As little shadows as possible)
4) The light should reflect from the pattern to the camera (as much as possible).
Thus, objects made up of metal are not suggested.
5) There should be minimal blur. If focus/depth-of-field is an issue,
use a smaller plate or less tilts.
6) If you can't fix the camera, it shouldn't move at all.
For extrinsic calibration, moving the camera after calibrating will ruin the entire
calibration. Be careful about this as even small shakes or movements distort the accuracy.
7) Don't pre-process the calibration images.
Even if you do very smart things, you might perturb the accuracy.
8) Ideally, have 15-50 images fitting the above criterion.
This should give you 450-1500 points to optimize in total.
9) If you are printing the pattern on a piece of paper, then make sure that
contrast between black-white squares is significant and no strange shape is appears on the paper.
10) Don't bend the calibration pattern. Make sure that it is as rigid (and planar) as possible
11) If you are printing on an A4 paper, make sure you have some white space around the pattern,
remember you will want to hold it somehow.
12) If you have distortion, use small circles, they are less affected from distortion.
If you have low resolution but less distortion use larger circles.
13) Start with low tilts and fronto-parallel images.
Move to the tilts and rotations in all axes, gradually, covering the entire visible space.
If you have bias in the space coverage or distribution of points,
then it means that some regions of the space will be better calibrated.
Keep this as equal as possible. Some variation can be tolerated.

I hope these help.


Implementation of halfband pass filter I need to interpolate a complex valued bandlimited periodic function using local interpolation.

I can have the signal sampled at any frequency I want over at equispace intervals.

I am aware that for bandlimited periodic functions the best interpolation method is global interpolation (FFT).

However, for this project I need to use local interpolation.
I managed to perform the interpolation using Lagrange polynomials of order 10 but I require an oversampling of 2 to achieve the desired results.
My goal is to reduce oversampling to any extend even if I need a longer interpolator.

I recently read as much as I could about FIR filters but somethings are not clear to me.

Since the interpolation factor is of the order of 2, I am planning to use a halfband pass interpolator (many books say this is the best since one gets a cheap interpolator due to the zeros in the impulse response).

My questions are:

How do I design a halfband pass filter? In Matlab there are nice functions in the dsp package that allow directly to construct such filters, unfortunately I have to do this implementation in C. I searched and I found that the remez() function from Octave to compute the coefficients is available in C and is the one I am planning to use.
Designing a low pass filter with the Parks-Mclellan algorithm is it by default a halfband pass filter (provided I set up correctly the frequency bands) or do I have to enforce the time domain constraints?
Regarding the linear phase of the FIR filter, what I understand is that my signal will be shifted by some samples, how do I compensate to remove the shifting? What comes to my mind is to apply the corresponding shifting on the opposite direction, is this necessary?
If you could provide an working example of an interpolation of a periodic function (e.g. sum of sines) using a FIR I would be very grateful.

The impulse response and frequency response of a length N zero-phase FIR lowpass filter can be designed using the Remez exchange algorithm.

N = 11; % filter length - must be odd

b = [0 0.1 0.2 0.5]*2; % band edges

M = [1 1 0 0 ]; % desired band values

h = remez(N-1,b,M); % Remez multiple exchange design

The impulse response h is returned in linear-phase form, so it must be left-shifted (N-1)/2 = 5 samples to make it zero-phase.
Answer taken from https://www.researchgate.net/post/How_can_I_design_a_null_phase_FIR_filter.

I implemented this solution some time ago and still works very well in my code.


What kind of filter would have a sawtooth frequency response? I want a complex FIR or IIR filter whose response is approximately a sawtooth in the frequency domain. How should this be formally described, and what formula or algorithm can I use to calculate taps for one?
Here is a sketch of the magnitude response I want (not analysis of an actual filter). I think the phase response doesn't matter.

A windowed design is easy. Starting with an ideal zero-phase frequency response:

$$\frac{\omega}{2\pi} + \frac{1}{2}$$

The inverse Fourier transform of that gives the ideal impulse response:

$$\int^\pi_{-\pi}\left(\frac{\omega}{2\pi}+\frac{1}{2}\right)e^{i\omega k}d\omega = \frac{\sin(\pi k)}{k} + i\left(\frac{\sin(\pi k)}{\pi k^2} - \frac{\cos(\pi k)}{k}\right).$$

Similar to the $\text{sinc}$ function, the real part $\frac{\sin(\pi k)}{k}$ has an illegal division by zero at $k=0$ where the limiting value $\pi$ should be used instead. For other integer $k$, the real part is zero. The imaginary part has similar problems at $k=0$, where its limit is 0. Drawing from Rick's answer, at integer $k$, $\sin(\pi k)$ is zero and $\cos(\pi k)$ is +1 or -1. This simplifies the ideal impulse response at integer $k$ to:

$$\begin{cases}\pi&\text{if }k=0,
\\i\frac{-(-1)^k}{k}&\text{otherwise.}\end{cases}$$

If the filter is implemented as a sum of a purely real and a purely imaginary filter, the purely real filter requires no work other than to multiply the input by a constant. The imaginary part of the impulse response is antisymmetrical, which may enable savings.

To continue with the windowed design, choose a window function and its length (equals the number of filter taps) and multiply the ideal impulse response by it.
That gives the impulse response of the filter. The window function should be centered at $k = 0$.

I plugged a frequency response of $\omega$ into the inverse discrete-time Fourier transform (inverse DTFT) equation and ended up with:

$$h[k] = \frac{sin(k\omega_c)}{k^2\pi}-\frac{(\omega_c)cos(k\omega_c)}{k\pi}$$

where $\omega_c$ is the cutoff freq measured in radians.  Index $k$ is ...-3,-2,-1,0,1,2,3....  When $k = 0$ set $h[0] = 0$. (You may want to multiply $h[k]$ by a window sequence to reduce the frequency magnitude response ripples.) If we set $\omega_c = \pi$ we end up with:

$$h[k] = \frac{-(-1)^k}{k}$$

Try out the above $h[k]$'s.  Who knows, they might work for you. (The $h[k]$ equations are supplied "as is." No warranty of any kind, either expressed or implied.)


Meaning of Real and Imaginary part of Fourier Transform of a signal Say $f$ is a signal of time $t$, $F$ its Fourier transform of the variable $v$.

It is known that in polar coordinate, $|F(v)|$ tells us how much the frequency $v$ is present over the signal, and $Arg(F(v))$ tells us how much the contribution of this frequency is phase-shifted.

What information do its real and imaginary part tell us ?

Or if I reformulate my question : can we give an interpretation of the Fourier transform in Cartesian coordinate like we can do in polar coordinate ?

The real and imaginary parts of the Fourier transform of a signal $x(t)$ are the Fourier transforms of the signal's even and odd parts, respectively:

$$X_R(\omega)=\frac12[X(\omega)+X^*(\omega)]\Longleftrightarrow\frac12[x(t)+x^*(-t)]=x_e(t)\\
X_I(\omega)=\frac{1}{2j}[X(\omega)-X^*(\omega)]\Longleftrightarrow\frac{1}{2j}[x(t)-x^*(-t)]=-j\cdot x_o(t)
$$

where $X_R(\omega)$ and $X_I(\omega)$ are the real and imaginary parts of $X(\omega)$, and $x_e(t)$ and $x_o(t)$ are the even and odd parts of $x(t)$, respectively.

Fourier transform of a system is its transfer function and gives the multiplication factor when $e^{j\omega t}$ is  the input.  $\omega$ is frequency.  If you consider the input as current, the transfer function or Fourier transform as impedance then the output is potential. If Fourier transform is impedance, then the real part of FT is resistive part of the impedance and imaginary part is the reactive part of the impedance.


Opposite of wavelet transform? Wavelet transform gives good time resolution for high-frequency events and good frequency resolution for low-frequency events.

=> I want to have complete opposite of wavelet transform, where I get good time resolution for low frequency events? Is there any known transform or specific window to get it?

wavelet transform also serves to you what you want to do, all you have to do is to first apply a low-pass your signal so as to keep the frequency range you want to scan only, another approach is to transform Gabor where you have to define the size and shape of the analysis window, but I would recommend better  using wavelet transform applying a lowpass filter your signal first.


How to understand the type of noise Suppose that I am receiving RF signals contaminated with noise through antenna.
This signal is digitally sampled.
How can I understand the type of noise(Gaussian,Uniform etc).
Any algorithm is there to test the type of noise.

There are several different ways to look at noise:

The probability density function (PDF) of a noise process, $x(t)$ say. It is sometimes written $p_x(x)$.  This tells you how frequently the the variable $x(t)$ takes on specific values.

One way to think of a PDF is as a (normalized) histogram of the possible values that $x(t)$ can take.  The PDF is normalized because $x(t)$ must take on some value between $-\infty$ and $+\infty$, so the integral
$$
\int_{-\infty}^{+\infty} p_x(x) dx = 1.
$$

For example, the figure below shows the PDF (in blue) and the histogram (in black) for a normal (a.k.a. Gaussian) random variable and for a uniform random variable.

Another way to look at noise is whether there is a relationship between the value of $x(t_1)$ and the value $x(t_2)$ (where $t_1 \not= t_2$). If there is no relationship, then the noise is said to be "white".  If there is a relationship, the noise is said to be "colored".

The way this is usually illustrated is using autocorrelation (autocovariance) of the signal.

Below are two plots of the autocorrelation of normally distributed noise.  In the first plot, the noise is white.  In the second plot, there is a high level of correlation between subsequent samples.

scilab CODE BELOE

// 15426

X_gauss = rand(1,1000,'normal');

X_uniform = rand(1,1000,'uniform');

figure(1)
clf
subplot(211)
histplot(100,X_gauss)
vals = [-3:0.1:3];
plot(vals,exp(-vals.^2/2)/sqrt(2*%pi))
subplot(212)
histplot(100,X_uniform)
plot([0 1], [1 1])

figure(2)

X_white = X_gauss;
X_colored = filter(ones(1,50),10,X_white);

figure(2)
clf
subplot(211)
plot(xcorr(X_white))
mtlb_axis([970 1030 -1000 1000])
subplot(212)
plot(xcorr(X_colored))
mtlb_axis([970 1030 -1000 1000])

The way I'm familiar with is:

1- Send the RF Image multiple times. Say 100 times.
2- Compare by "subtraction" in case of an Image for example .. between the received image/ and transmitted one.
3- The subtraction will give you the noisy image.
4- Plot a histogram statistics for the noisy image, for every received image. You'll then be able to roughly model the noise based on the experimental results.

The most familiar ones are, Guassian, or even Uniform.

if it's not of these familiar ones. Run a regression code to represent this model as an equation. By that means you can take the inverse of it to compensate for the noise.


Calculation of an autocorrelation function A sample of a random process is given as:

$$
x(t) = A\cos(2\pi f_0t) + Bw(t)
$$

where $w(t)$ is a white noise process with $0$ mean and a power spectral density of $\frac{N_0}{2}$, and $f_0$, $A$ and $B$ are constants.
Find the auto-correlation function.

Here's my attempt at a solution:

Let $a = 2\pi f_0t$, and $b = 2\pi f_0(t+\tau)$

\begin{align}
\text{Autocorrelation of } x(t) & = E\left\{x(t)x(t + \tau)\right\}\\
& = E\left\{\left(A\cos(a) + Bw(t)\right)\left(A\cos(b) + Bw(t+\tau)\right)\right\}\\
& = E\{A^2\cos(a)\cos(b) + AB\cos(a)w(t+\tau) + AB\cos(b)(wt)\\&\quad + B^2w(t)w(t+\tau)\}\\
& = E\left\{A^2\cos(a)\cos(b)\right\} + E\left\{AB\cos(a)w(t+\tau)\right\} + E\left\{AB\cos(b)(wt)\right\}\\&\quad + E\left\{B^2w(t)w(t+\tau)\right\}\\
& = E\left\{A^2\cos(a)\cos(b)\right\} + E\left\{B^2w(t)w(t+\tau)\right\}\\
& = E\left\{A^2\cos(a)\cos(b)\right\} + B^2\left(R_w(\tau)\right)\\
& = E\left\{A^2\cos(a)\cos(b)\right\} + B^2\left(\frac{N_0}{2}\right)(\delta(\tau))\\
\end{align}

The expectation terms with the noise in them all equal $0$ (the last is just the auto correlation of white noise ... hence the simplification above. Using trigonometric identities:
$$
\cos(a)\cos(b) = \frac 12\left[\cos(a + b) + \cos(a - b)\right]
$$

we have:

\begin{align}
\text{Autocorrelation of } x(t) & = E\left\{A^2\cos(a)\cos(b)\right\} + B^2\left(\frac{N_0}{2}\right)(\delta(\tau))\\
& = E\left\{\left(A^2\right)\frac 12\left[\cos(a+b)+\cos(a-b)\right]\right\} + B^2\left(\frac{N_0}{2}\right)(\delta(\tau))\\
& = \left(\frac{A^2}{2}\right)\left[E\{\cos(a+b)\} + E\{\cos(a-b)\}\right] + B^2\left(\frac{N_0}{2}\right)(\delta(\tau))\\
\end{align}

We're dealing with constant terms, so expectation term goes away and subbing in our initial conditions we get:
$$
\frac {A^2}2 \left[\cos(2\pi f_o(2t + \tau) + \cos(2\pi f_o\tau)\right] + B^2\left(\frac{N_0}{2}\right)(\delta(\tau))
$$

For some reason I can't help but feel I did something incorrectly calculating that autocorrelation ... it's supposed to be a function of $\tau$, but has a $t$ is in there ... I would very much appreciate it if someone could point me in the right direction, or explain what I messed up. I don't know whether it matters, but in this class we're dealing with only wide sense stationary processes.

I guess you've done almost everything right, but have a problem at the calculation of the expectation value regarding $t$. You should calculate the expectation value of the cosine function. Sadly, it does not simply "go away" as you wrote.

Have a look at the Wikipedia page. There you can find another, more explicit, formula for the auto-correlation function of a function $f(t)$:

$R_{\textrm{ff}}(\tau) = \int\limits_{-\infty}^\infty f(t+\tau)\bar{f}(t) \, \textrm{d}t$ .

(Note, that compared to the Wikipedia page, I have taken the liberty to use the variable $t$ in the integration instead of $u$, which would be the mathematically more accurate version.)

As you can see from this equation, you "integrate away" the dependency on t, and indeed you should be left with a function that is independent of $t$.

Note that there is also an version that does not go to infinite times, but is constrained to a period $T$. Maybe this version is more appropriate in your case. However, the same holds for this version: $t$ is integrated away and should not be a variable in the resulting formula.


Log Distance-Normal Shadow Model and Friis Since I'm founding several but different results about this model, I would ask you if all it is OK with this formulas:

1) Friis Model --> $Pr_0 = Pt + Gt + Gr + 20\log_{10}\frac{\lambda}{4\pi d}$

2) Log Normal Shadow --> $Pr_L = Pr_0 - 10\eta\log_{10} d + \chi_{\sigma}$

where:

$Pr_0$ is the power received in the Friis Model as depicted above
$Pr_L$ is the power received in the Log Normal Shadow Model as depicted above
$Pt$ is the power transmission
$Gt$ and $Gr$ are the respective transmitter and receiver gains
$\lambda = c/f$, where $c$ is the light speed and $f$ is the frequency
$\eta$ is the Path Loss exponent
$\chi_{\sigma}$ is a Gaussian zero mean random variable with standard deviation $\sigma$.

I am a bit concerned about the signs and the $\log$ operations.

I verified all the steps and they are all correct.


Adaptive volume normalization on a stream of PCM data I am writing a program that receives a continual stream of PCM data from various sources. Some sources have different volumes than others. The PCM data is always a live continuous stream of 16 not signed data single channel.

How can I dynamically make sure that all samples are not too low in volume or high in volume to effectively normalize its output realtime?

What you want is called automatic gain control or dynamic range compression. It can be implemented in a dynamics processor that takes the input audio, uses an envelope follower to analyze the volume of the input as it changes over time and modifies the audio so that the envelope of its output corresponds to the input envelope put through a non-linear function such as:

Adobe Audition Dynamics Processing

I don't know what is the best practice in how the envelope should be modified to match the target envelope, but usually there are separate parameters for attack and decay that control how fast the applied gain should be increased or decreased, respectively, to make the resultant envelope approach the target envelope. You could try the potentially unstable:

$$\text{gain}(t) = \text{target_envelope}(t)/\text{input_envelope}(t)$$

where $t$ is time. The resultant $\text{gain}(t)$ should be smoothed (lowpass filtered) before applying it, to get the wanted attack and decay behavior. The envelope detector can have additional controls for changing the (usually root mean square or peak) analysis window, and a lookahead time setting that enables the dynamics processor to respond to volume changes before they take place, by delaying the input audio after envelope detection. The root mean square calculation can be done by calculating the square of the input audio, lowpass filtering that and taking the square root of the result.

Adobe Audition Dynamics Processing attack, decay and lookahead controls

These parameters and the non-linear function should be tuned (by listening) for minimal audible artifacts such as pumping of the volume and low-frequency distortion.


What level of correlation coefficient suggests correlation? In Python (and Numpy), I recently did:-

a = np.random.randint(0,256,1000)

b = np.random.randint(0,256,1000)

and when I calculate the correlation coefficient $r$ using:-

np.corrcoef( a,b )

I get a value of 0.07944437.  Rounded, that's 0.08, 8% or about $1 \over 13$ th.

$r$ is said to typically range -1.0 thru 0.0 to +1.0.  I've never had a 0.0 and don't think that it's possible as random correlations can always be found between any two series. It's like asking "How small is small?" 0.08 doesn't seem that small. Yet clearly $a$ is not correlated to $b$ as both originate from the good quality Mersenne Twister deep inside Python.

So if $r = 0$ is empirically impossible, when does it actually become significant and suggest a real correlation?

Note 1.  There is of course the 68–95–99.7 rule. This covers a $3 \delta$ range of certainty. The atomic guys use $5 \delta$ to confirm a discovery. Do  signal guys do anything similar, and if so, how would they apply it?

Note 2. With the greatest of respect, How do I know quantitatively if the correlation of two time series is significant? seems ambiguous.

One way to look at it is that $r^2$ (yes, I do mean $r^2$, not $r$ or $|r|$) is the fraction of the variance of the target (the series A) that is explained by the estimator (a linear function $\alpha + \beta B$ of the series B).  Note that if Series A had variance $\sigma^2$, then the variance of $(A - (\alpha + \beta B))^2$ (what B cannot explain about A) is $(1-r^2)\sigma^2$.  So, if $|r| < \frac{1}{\sqrt{2}} \approx 0.707\ldots$, less than half the variance has been explained, and so one viewpoint is that a correlation of less than $70\%$ or so suggests that the linear estimator being considered here is not all that great.


Parallel equalizer, working for boosts but not cuts I am trying to replicate the behaviour of a parallel equalizer.

The approach is to use a parallel sum of bandpass filters like so:

$Response = 1 + bandpass1 * (gain1 - 1) + bandpass2 * (gain2 - 1) + ...$

Boosts and cuts are separately grouped and run in series.

I have successfully got the model working for two boosts in parallel. But running two cuts in parallel produces an unexpected result.

The figure below illustrates the problem. The thick blue line illustrates two boosts in parallel (thin orange and red line show the individual boosts that are being combined). The thick purple line illustrates two cuts in parallel (thin green and light-blue line show individual cuts that have been combined). As can be seen the individual cuts are a mirror image of the individual boosts but the parallel sum doesn't work for cuts.

Does anyone know what I'm missing so that that the thick purple line is a mirror image of the think blue line?

Here as my matlab code:

% boost 1
f1 = 176;
g1 = db2gain(12);
q1 = 1;

% boost 2
f2 = 672;
g2 = db2gain(12);
q2 = 1;

% cut 1
f3 = 176;
q3 = 1;
g3 = db2gain(-12);

% cut 2
f4 = 672;
g4 = db2gain(-12);
q4 = 1;

% bandpasses for boosts
[b1, a1] = bandpass(f1, proportionalQ(q1, g1));
[b2, a2] = bandpass(f2, proportionalQ(q2, g2));

% bandpasses for cuts
[b3, a3] = bandpass(f3, proportionalQ(q3, g3));
[b4, a4] = bandpass(f4, proportionalQ(q4, g4));

f = logspace(log10(20), log10(20000), 10000);

r1 = response(b1, a1, f); % bandpass response for boost 1
r2 = response(b2, a2, f); % bandpass response for boost 2
r3 = response(b3, a3, f); % bandpass response for cut 1
r4 = response(b4, a4, f); % bandpass response for cut 2

% run the analog bandpass filters in parallel as two boosts (OK)...
p1 = 1 + r1 .* (g1 - 1) + r2 .* (g2 - 1);

% run the analog bandpass filters in parallel as two cuts (NOT OK)...
p2 = 1 + r3 .* (g3 - 1) + r4 .* (g4 - 1);

figure;

semilogx(f, response2db(p1), 'LineWidth', 2); % parallel boosts (OK)

hold on;

semilogx(f, response2db(1 + r1 .* (g1 - 1))); % individual boost plot 1
semilogx(f, response2db(1 + r2 .* (g2 - 1))); % individual boost plot 2

semilogx(f, response2db(p2), 'LineWidth', 2); % parallel cuts (NOT OK)

semilogx(f, response2db(1 + r3 .* (g3 - 1))); % individual cut plot 1
semilogx(f, response2db(1 + r4 .* (g4 - 1))); % individual cut plot 2

xlabel('Frequency');
ylabel('Magnitude (dB)');

axis([20, 8000, -16, 16]);
grid;

hold off;

bandpass.m

function [b, a] = bandpass(f, q)
    b = [0, f / q, 0];
    a = [f * f, f / q, 1];
end

response2db.m

function [db] = response2db(r)
   db = 20 .* log10(abs(r));
end

db2gain.m

function [g] = db2gain(db)
    g = 10 .^ (db ./ 20);
end

proportionalQ.m

function [q] = proportionalQ(q, g)
    q = q .* sqrt(g);
end

response.m

function [r] = response(b, a, f)
    s = 1i .* f;

    r = (b(1) + s .* b(2) + s .* s .* b(3)) ...
        ./ (a(1) + s .* a(2) + s .* s .* a(3));
end

Cuts can simply be derived from boosts like so: p2 = 1 / p1

This is intuitive as this is how you would invert a bandshelf.

However, I'm not sure it leads to an elegant equation as inverting the sum of filters and recovering a rational filter expression might be messy, so there's potential a better way.

I can't run your code since one function is missing.

In general, parallel EQs are a not a great idea. Since you add the responses, the exact phase responses of either boost or cut need to be carefully controlled, otherwise you see cancellation (as you probably do here).

It's much easier to parametric EQs in series instead. If the gain is 0, there is no interference in they add simply in dB.

This being said : this line looks strange to me

p2 = 1 + r3 .* (g3 - 1) + r4 .* (g4 - 1);

I assume g3 is supposed to be the amount of cut. You get no cut for g3 = 0dB but it's certainly NOT -12dB for g3 = -12dB.  (gr-3) comes out to be -.75  and gives you a nasty phase flip as well


Efficient real-time pitch detection algorithm What is the most efficient pitch detection algorithms in terms of computation? I need something to implement on a DSP.
The typical signal in my situation would be a mixture of low-level wide-band noise plus a series of strong harmonics below 1000 Hz (think of the noise produced by industrial machinery or by car engines). The system needs to sample ("know") the current fundamental frequency of the signal every second or so.

YINS fft is the best method that worked for me in the past. I have my open source implementation of it written in awesome D language: link . I compared the results with Praat they match.


FFT: Removed padded zeroes? I have an application where I do this:

DFT->Filter->IDFT on a range

For computational performance I'm zero-padding my FFT to a power of two, but when I take the IDFT I want to get out an interpolated image instead of the zero-padded image. A minor detail, is that the IDFT is performed over only a certain frequency window.

What can I do in the Fourier domain to stretch my signal so that the IDFT returns the image without the zeroes I padded it with?

To make it clear, you padded your image in what domain ?

Time domain >> No problem, IFFT and then truncate will bring you the original image. It's still intact.
Freq domain >> Padding in the Freq domain means image resizing. You cant get your original image with 100% anymore. If you really want the original size one, the only thing you can do is to resize it back with some resampling methods.

First of all, common FFT libraries (which software are you using?) have FFT implementations for a lot more lengths than powers of two -- the FFTW, which is by far the most widely used in general purpose computing (and is in use by things like Matlab, FFMPEG and so forth) for example supports any length that only has prime factors smaller than 13 or so; so unless you're doing an FFT of millions of points in each direction, you don't have to pad for performance reasons.

I assume you pad before transforming:

signal -> padding -> DFT -> Filtering -> IDFT -> ?

Now, you can just truncate your result, and will get the right thing -- the IDFT works cyclic, and thus you won't be throwing away information.

Also note that if "Filter" in your original description is a FIR filter, you are convolving your signal in frequency domain, just to transform it back --- just multiply it with the fourier transform of your Filter, and you'd have done the same in the original domain of your signal, as basic Fourier theory allows you to do.


making pink noise (1/f) using list of frequencies making pink noise (1/f) using list of frequencies

I would like to see what type of noise I would get if I used just the frequency in my voice.  I created a matlab array using fft to get the
[frequency,amplitude,phase] to reproduce my vocal signal.

I would like to take this file/data and use it to create pink noise (1/f).  Of course when I use 1/f for the frequency the numbers become very small does anyone have any ideas how to use my own vocal frequencies that I get from doing a fft in matlab from my voice to create pink noise (1/f).

Thanks

Pink noise, as you know, has contributing amplitudes inversely proportional to each corresponding frequency. It appears more sinusoidal than white noise (constant power spectrum) because the low frequency components have large amplitudes. The high frequency components do not contribute that much to the overall noise amplitude.

To recreate pink noise in MATLAB, I would use the following technique:

1.) Generate a square wave (ones(n) will create a square wave of width n).

2.) Take the Fourier Transform of the signal (fft(x)). This will give you the sinc function. The DC component will be 1.

3.) In a for loop, add noise to the amplitude of each frequency of the FFT sinc function. The range of values of the noise will be dependent on 1/f.

4.) Take the inverse Fourier Transform of the noisy FFT. Then subtract 1 to yield a pink noise vector of length n centered around $0$.


MATLAB Hamming window I have been trying to achieve a Hamming window length of 23 using the equation.
I tried several coding, one of the is the one below:
M=23
w = .54 + .46*cos(pi*(-M:M-1)'/M)
plot(w)

If the length of the window should be $23$, $M$ must equal $(23-1)/2=11$:

M = 11;
m = -M:M;
w = .54 + .46*cos(pi*m/M);
plot(m,w)

There is a default function in MatLab which can generate this for you, called hamming. In case you do not have access to this function, because you do not have the signal processing toolbox, you could try:

M = 23;
w = 0.54 - 0.46 * cos(2 * pi * linspace(0, 1, M));


How to take acoustic measurements of a clasroom I'm an electrical engineering major and I have a question regarding a project I'm working on.

I'm supposed to design 3 band-stop filters on MATLAB for echo cancellation/suppression.

My problem is not filter designing, but with taking acoustic measurements of a classroom since I have little or no information in that area.

The supervisor suggested taking coefficients using impulse response invariant! I could not find something solid or a method while searching online.

  So my question is, what is the simplest way of taking acoustic
  measurements of a normal size room/classroom? How to find the
  coefficients through impulse response invariant method?

All suggestions are appreciated, current tools I have (dynamic mic, speakers, audio interface (Focusrite Scarlet 2i2), RoomEQ wizard, MATLAB/Simulink, FPGA DE2 board with cyclone)

I tried using RoomEQ wizard and followed their tutorial but did not work properly or maybe I`m not getting the whole picture.

It doesn't matter where your microphone is in the room when you take the measurement. What you need to do is find where the peaks are in the cross correlation (each peak corresponds to an echo at a certain lag/delay), figure out what the corresponding inverse filter needs to be to cancel out the sound.

If you're trying to do echo cancellation in MATLAB and have the signal processing toolbox, take a look at the example in xcorr by typing:

  doc xcorr

for a simple example.

Edit 1:To find the inverse filter required, you need to calculate the transfer function between the source and the mic. (edit: sentence in here was misleading) Once your filters are set up, as long as nothing is changing you should be able to cancel the echoes produced by the room.

Edit 2: To get the room impulse response you typically use a long LFM chirp (some people pop a balloon) and a microphone, then calculate the transfer function between the input signal and the recorded signal. Here's an example paper of what you can do: http://www.artalabs.hr/papers/im-aaaa2003.pdf

Edit 3:

  You do realise that reverberation (or echoes as you described) are dependent on location and even time?

you said you are making an "impulse response invariant" filter. So as long as neither your source or microphone are moving, the IIR filter created from your Room Impulse Response measurement will cancel out your echoes...just make sure the mic and speaker stay in the same spot.


How to perform iris normalization? I have done with iris localization. Now I need to perform iris normalization, but I don't have any idea about it. How to implement iris normalization please?

Please check the following for Iris Normalization:

  Nithyanandam. S, Gayathri K.S., Priyadarshini P.L.K. - A New IRIS Normalization Process For The Recognition System With Cryptographic Techniques

  Arvacheh E. - A Study of Segmentation and Normalization for Iris
  Recognition
  Systems

you can also this two papers for iris and normalization:

S Zeriouh, M Boutahri, SE Yamani, A Roukhe. "Contribution of Artificial Neural Networks to the Identification and Detection of Targets Concerning Mobility on Remote Sensing Images." International Journal of Computer Science and Software Engineering 4, no. 5 (2015): 130-136.

How to implement iris normalization?

R Safdari, Z Azadmanjir. "Wireless Networks in Hospitals; a Framework of Use and Implementation for Nursing Care." International Journal of Computer Networks and Communications Security 2, no. 6 (2015): 188-196.


A Practical Method for Correcting Phase Errors in a Transceiver System Comprising of OFDM + DUC and DDC I am back with a question related to my OFDM Transceiver System that I am designing where I use DUC and DDC. I could get rid of the high BER problem. I realized that the problem was due to the phase error occuring between the transmitter and the receiver side as I saw the constellation diagram at the receiver rotating constantly. So to get rid of this, I obtained the phase difference between the transmitter and the receiver and I rotated the complex signal with the obtained phase error. However, I feel this is only an ideal method to get rid of the phase error. I read about the phase error compensation techniques but those are mostly when there is a noise in the channel. But in my case, I still haven't added any noise in the channel. I suspect the phase error is either because of the DUC and the DDC, or because of the pulse shaping filters that I used after OFDM modulator.

So I was wondering if anyone could tell me why is this phase error occuring even without a noisy channel, and what can I do to mitigate the phase error with some estimation.

I have attached the screenshot of my system along with the waveforms obtained. Below are the specs -

Random Integer Generator Specs

Set Size - 2

Initial Seed - 37

Sample Time - 1/1e6

Samples per frame - 20

OFDM Modulator Specs

FFT length - 32

Number of guard bands - [6;5]

Insert DC Null - Checked

Cyclic Prefix Length - 16

Number of OFDM symbols - 1

Number of transmitted antennas - 1

Raised Cosine Transmit Filter Specs

Rolloff factor - 0.2

Filter span in symbols - 8

Output samplers per symbol - 4

Linear amplitude filter gain - 1

DIgital Up- Converter Specs

Interpolation factor - [2 8]

Minimum order filter design - Yes

Two-sided bandwitdth - 2.4 e6

Passband ripple - 0.05 dB

Stopband ripple - 90 dB

Center Frequency  - 30e6

Input sample rate - 9.6e6

Digital Down-Converter Specs

Decimation factor - [8 2]

Minimum order filter design - Yes

Two-sided bandwitdth - 2.4 e6

Passband ripple - 0.05 dB

Stopband ripple - 90 dB

Center Frequency  - 30e6

Input sample rate - 153.6e6

Thank you in advance!

I don't think it would be reasonable to expect the phase to be exactly the same. The raised cosine filter for example has a linear phase response. Also don't forget that the auto-correlation for the cyclic prefix detection might find the frame starting at any of the 16 samples of the cyclic prefix, inducing an undetermined phase shift.

Do you have pilot sub-carriers? One simple way is to analyse their phase (which you know a-priori were transmitted with 0 (or some other known) phase) and rotate the sub-carriers in between the pilots (linearly interpolating between the phase of the pilots).

I don't think OFDM uses a raised cosine filter, that breaks orthogonality.


MU-MIMO Sensitivity to Line of Sight (LOS) Why do researchers talk about MU-MIMO being "less sensitive to line of sight" / "LOS Okay"? For example in this video, time=[00:40,00:42]

Similarly, "User separation based on spatial channel properties is
particularly  difficult  in  situations  where  the  users  are  located
close to each other and experience LOS propagation conditions
to the BS antenna array"

My perception is that LOS is always preferred to the opposite, multipath propagation.

Can someone please explain this? Thank you.

There is no mention of light-of-sight (LOS) on the time you mentioned. There is a mention of point-to-point, however. Point-to-point channels could be LOS or multipath.

MU-MIMO is better than point-to-point MIMO because users are usually well-separated in space. MU-MIMO is particularly better when the channels are correlated at the base station (BS). The reason for that is that the spatial separation of the users make the channel correlation matrices from different users to be different. In the ideal case where these correlation matrices eigenspaces are orthogonal, the data from different users can be perfectly separated without any interference. In point-to-point MIMO, on the other hand, the data from different antennas from the multiple-antenna transmitter will have the same correlation matrix at the BS, and thus there will be more interference from other data streams.

Read section 2.4 from the book Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency.

EDIT: In the paper, they are talking about separating users that have LOS channels with the BS in massive MIMO, that are difficult to separate in traditional MIMO.

For MIMO to have any advantage, you need uncorrelated paths between the multiple antennas.

If all the users have LOS, then their path is just a complex number with no way to become uncorrelated; that makes the channel matrices non-separable using MIMO methods.

In other words: MIMO uses the fact that different antennas see different channel realizations. When effectively all channels are the same, good channel, then you don't get to use these differences.


EEG signal processing I have an EEG data set downloaded from physionet. I want to play that data into a brain monitoring device like for example BIS monitor or Narcotrend or any other similar brain monitoring device.I use Matlab GUI to visualize the data. Can I send it from a Matlab GUI to any one of the brain monitoring devices with a suitable hardware  in between.Any ideas please I am new to EEG signals.

I'll assume that you want to do online monitoring.

You can always store the coming data from the GUI in a vector in MATLABA.
Then you can run simple code to plot them.

This also depends on the type of measurements you have. If the set is capable on detecting power ranges as well you should have the following readings:

High Alpha, Low Alpha, High Beta, Low Beta, High Gamma, Low Gamma, Delta, Theta, Attention, and Meditation.


Equalize Speaker using Frequency Response Curve I have a speaker for which I know the frequency response curve (-5dB at 600Hz, -3dB at 700Hz, ... +5dB at 10kHz). This response curve is not as flat as I would like, and so I want to write a script which takes an audio file as input, and then transforms it such that when this new audio is played out through the speaker, it sounds like what the original file would sound like if it were played through a speaker with an actually flat response curve. (Basically, I want to flatten the frequency response of this speaker in software.)

What would be the best way of going about this? My first instinct is to Fourier transform the audio signal, reweight each Fourier coefficient according to the inverse of the amplitude at the corresponding frequency in the frequency response curve, and then inverse Fourier transform the result. However I imagine that such a direct approach would result in problems such as ringing artifacts, and I know that there are ways of performing similar tasks through the use of specially designed filters.

Oops I misread your question and thought you have a microphone recording of the speaker output with given audio, and answered that. Well here is anyway, as some parts still apply:

Like others have suggested in comments, designing a filter is the way to go. I would recommend a windowed FIR filter design rather than Parks-McClellan. From MATLAB documentation:

  firpm exhibits discontinuities at the head and tail of its impulse
  response due to this equiripple nature.

Our ears analyze sound with kind of wavelets so they may detect those discontinuities before and after impulsive sounds. With windowed filter design however you can choose a window that smoothly starts from and ends at value 0, avoiding creating such discontinuities.

Doing division (inverse of the amplitude) is a bit dangerous, because if you divide something by a number close to zero you get a huge value. Discrete Fourier transform results often contain values close to zero. If you have a long audio file then you can split the audio into shorter pieces with say 50 % overlap, multiply each with a window function, and calculate the mean squared magnitude of their discrete Fourier transforms (discarding phase). The square root of this is the root mean square spectral power, which is unlikely to be zero anywhere. Do this for both the speaker and raw audio and bin-wise divide the root-mean-square spectral power of the latter by the first. Still, you may get very large values if the speaker is unable to reproduce very low or very high frequencies. You can set a (soft) limit to how much correction you do at most, and you can even softly ramp the correction to nothing at frequencies you know the speaker simply cannot reproduce. Assuming zero phase for all bins, you can directly inverse discrete Fourier transform the division result to get an unwindowed zero-phase impulse response. You would then multiply the unwindowed impulse response with a window function which can be much shorter than the transform result, giving you a practical, short filter impulse response.

One more trick is to convert the filter to its minimum-phase counterpart, which will completely eliminate pre-ringing (taking place before the main peak in the impulse response). There are specialized methods for doing this that preserve the magnitude frequency response. I would not worry about post-ringing because real-world acoustics will give you plenty of post-ringing and echos anyhow and we are used to masking them.

Couple software you could use.

REW - http://www.roomeqwizard.com/

REW inputs the speaker FR data (measured or entered) and lets build the correction filters/data to be used with EQ software (EqualizerAPO as for an example). When you have the correction data in hands there are many methods you can use to get the correction embedded into your audio file (if this is your aim to do?).

EqualizerAPO - https://sourceforge.net/projects/equalizerapo/

I usually make the inverse of the FR curve by loading the response curve plot (as certain size image file) as a background image to my own fPEQGUI-10 (EqualizerAPO GUI software) then adjust filters (10 various type filters available) output to match the imported FR curve and then finally inverse the filters to get the correction filters/coefficients out (few output formats supported). Benchmark utility (bundled with EqualizerAPO) can be used for to generate the audio file (or IR file) if needed.


Z transform stability What is the causality & stability status for three cases shown (aso in attached photo) ?

$$H(z) = \frac{z(z-1)}{(z+1)(z+\frac{1}{3})} $$

for three possible regions of convergence as:

a-) |z| > 1

b-) |z| < 1/3

c-) 1/3 < |z| < 1

Am i ok as below?

a) Causal & marginally stable

b) anti-causal & unstable

c) non-causal & marginally stable

Yes, all answers given by you are fine.

Causal: If ROC is outside the outermost pole

Left sided: If ROC is inside the innermost pole

Stable: If ROC contains the unit circle (marginally stable if it touches unit circle)

I will only give you hints
1. Casual if |Z|>|a|
2. Stable if Roc contains unit circle
So non causal if |Z|<|a|, unstable if Roc don't contain unit circle & marginally stable if poles are on unit circle.


Do symmetric discrete signals have zero phase? I generated a Hanning window having an even symmetry:

where for even-sampled case we either take "left" and "right" to include or exclude the center sample. I was surprised to find its DFT has non-zero imaginary coefficients, i.e. non-zero phase, in both cases. I've also tried an odd symmetry case and ended up with non-zero real coefficients.
Shouldn't even- & odd-symmetric discrete signals have zero-phase DFT?

Assuming you are referring to the standard definitions of odd and even symmetry, i.e.
\begin{align}
\text{even: }&& x[-n] &= x[n] \\
\text{odd: }&&x[-n] &= -x[n]
\end{align}
then it still holds that the DFT of an even sequence is real (i.e. the phase is 0 or $\pi$) and the DFT of an odd sequence is imaginary (i.e. the phase is $\pi/2$ or $-\pi/2$).
This is independent of whether the the length of the sequence is even or odd as long as the sequence meets the definition of symmetry. This implies that for odd sequences $x[0]$ must be $0$ and for odd sequences of even length $x[N/2]$ must be zero too, where $N$ is the DFT length.
Both sequences that you show in your original question are not "symmetric" in the mathematical sense. They are only symmetric when properly time shifted and time shifting is clearly a non zero-phase operation.
The following Matlab script illustrates that this is true within the numerical accuracy of double precision math:
%% script to confirm zero phase DFT of even/odd symmetric sequences

N = 128; % FFT length

%% Even length sequences

% create a random sequence of length N/2+1
x = randn(N/2+1, 1);

% create even version
xEven = [x; x(end-1:-1:2)];

%  perform DFT and check zero phase
fx = fft(xEven);
zeroPhaseError = 10 * log10(sum(imag(fx).^2) / sum(abs(fx).^2));
fprintf('Even symmetry, Even length: zero phase error = %6.2fdB\n', zeroPhaseError);

%% Odd length sequences

% create odd symmetry version. This requires the samples at n=0 and n=N/2 to be
% zero
xOdd = [0; x(2:end-1); 0 ; -x(end-1:-1:2)];
%  check zero phase
fx = fft(xOdd);
zeroPhaseError = 10 * log10(sum(real(fx).^2) / sum(abs(fx).^2));
fprintf('Odd Symmetry, Even Length: imaginary phase error = %6.2fdB\n', zeroPhaseError);

%% now try an odd length sequence
% create a random sequence of length N/2
x = randn(N/2, 1);

% create even symmetry version
xEven = [x; x(end:-1:2)];

%  check zero phase
fx = fft(xEven);
zeroPhaseError = 10 * log10(sum(imag(fx).^2) / sum(abs(fx).^2));
fprintf('Even symetry, odd length: zero phase error = %6.2fdB\n', zeroPhaseError);

% create odd symmetry version. This requires the samples at n=0 to be zero.
% There is no sample at n=N/2 since that's not an integer.
xOdd = [0; x(2:end); -x(end:-1:2)];
%  check zero phase
fx = fft(xOdd);
zeroPhaseError = 10 * log10(sum(real(fx).^2) / sum(abs(fx).^2));
fprintf('Odd symmetry, odd length: imaginary phase error = %6.2fdB\n', zeroPhaseError);

TL;DR No, they don't. But in a specific sense, their continuous-time interpolation will be zero-phase if symmetric about $t=n=0$ (i.e. in periodic extension). In discrete case, however, symmetric will have approximately zero phase - the more samples the better.
All code. References: Wiki, MATLAB, scipy.

To see why no symmetric discrete signal can be zero-phase, begin by considering a fewer-sample random signal; the base sequence throughout the answer will be:
$$
s=[.4, -.4, -.2, .1, -.9, .2],\ N=6
$$
First, symmetric even and odd, formed by concatenating reflections as shown in titles, and their DFTs:

Now; how do we concatenate $s$ such that it is zero-phase? For best results, try to code this yourself. Answers:

 below

Not very symmetric. But why does it work?
Rooting the DFT are its complex sinusoidal bases, where the answer lies; for the even, left-right symmetry case, the real part is nonzero - i.e. cosines (recall what DFT coefficients mean). Consider $k=1$, even and odd (of same lengths as above waveforms):

Observe closely. Where is symmetry in either? And in both? For all-cosines to add up to something that is symmetric about some point, they themselves must be symmetric about that point. Where is that point?
Nowhere. There is no complete symmetry; this is an emergent property of the fact that the bases span one full period, i.e. we do not include the right-most $1$ in the cosine, as that overlaps into the next period. Only if we did would we have symmetry.
Then how do we get zero phase? Well, while there is no "clean" symmetry, there is in fact a point of symmetry in a less straightforward sense. Re-examining above plot, try and see if you can attain symmetry by excluding one point. Hint:

 Drop left-most point, $n=0$

Answer:

It's exactly at $N / 2$. Note for odd case we drew a line where no sample lies for a visual, but it is consistent with the definition of odd-sampled even symmetry in the question. Further, note that this holds for every other basis up to Nyquist, not just $k=1$.
With this information in mind, the reader is encouraged to revisit zero-phase plots above. Now to restate our findings:

Even-sampled zero phase: attained by symmetry about $N / 2$, ignoring sample $n=0$. This amounts to two "don't care" points, one at $n=N/2$ and other at $n=0$.
Odd-sampled zero phase: same, but only one don't-care. Visually this is a flatline as opposed to a spike in even case.

Now we're finished. All of above can be repeated for odd-symmetry case, i.e. pure-imaginary coefficients, or sine basis; the point of symmetry's the same for the same reason it's there for cosine. (However, odd symmetry for a signal is trickier as now the vertical axis is involved, so extra steps are due before concat, but this isn't a problem for zero-mean sines).

What about $n=0$?
Dissenters might point to this; it's readily invalidated by considering the even-length case - reproduced below, and fftshift-ed to center about $n=0$ on right:

It does work out for the odd-length case, however:

Now as to in what sense $n=0$ is always a point of symmetry; we must leap into the continuous. Recall again the reason why we have the kind of symmetry we have in the discrete case: single period. Now, in continuous, our bases still begin at $0$, but no longer end at $N-1$; instead, we have $t\in[0, N)$ - a semi-open interval, all $t$ from $0$ up to but excluding $N$.
So do we still have a point of symmetry at $N/2$? Yes - but what changes is the interval of symmetry, i.e. what counts as "left" + "right" has now shifted from spanning $[1, N-1]$ to spanning $(0, N)$:

(Lines are actually exactly at $0$ and $N$ but I shifted slightly for clarity). What does this change in regards to $n=0$ symmetry? Well, now we imagine we have a continuous-time signal spanning $[0, N)$, and it's easy to redraw earlier figures with this in mind and see that there is indeed a semi-open interval symmetry about $n=0$.
Above constitutes a fair definition of symmetry about $n=0$. Note, however, that this is not same as discrete, finite symmetry, and a discretely symmetric waveform in time domain about "center" ($N/2$) or generally about $n=0$ (i.e. accounting for odd case) will not have zero phase.
To test for zero phase, apply the exact criterion in the question not to $x$, but to $x[1:]$.


Biquad cascading frequency response I am implementing EQ lines in my filter. I need to do biquad with slope 24 db/oct, 36/oct and bigger so I would like to cascade several biquads. I am stuck with calculating frequency response of two biquads - now I am able to calculate magnitude and phase and draw one curve (peak filter).
Should I multiply this response by itself? Thank you in advance for the answers.

The application of a filter to a signal in the time domain is defined as convolution.  Convolution in the time domain is mathematically equivalent to multiplication in the frequency domain. So if we filter a signal, and then filter the output of the first filter, this is equivalent to multiplying the signal by both filters in the frequency domain. Therefor, this is equivalent to filtering with a third single filter whose frequency response is equal to the product of the first two filters’ frequency responses. Because frequency responses are complex numbers, the resulting magnitude is equal to the product of the two input magnitudes and the resulting phase is the sum of the input phases.  This can be easily extrapolated for more than two filters as well.

Thank you for the answer, but i need some clarification - i think i am missing something important to draw this plot.
I am doing steps below:

calculate in for loop freq response of standard biquad (peak) and get array of complex responses
get abs() from each response and it gives me magnitude for given frequency on plot
do 20 * log10(magn) to get values in dB and i am getting for example peak line with min in -9dB
Phase i am calculating for each magnitude in first biquad using atan2(imag(freqRes)/real(freqRes)).

Then i am trying to calculate 4th order so cascading the biquad:

each complex response from previous biquad i am multiplying by itself (freqRes)^2
an multiply it by exp(2jphase) so (freqRes)^2 * (cos(2jphase) + sin(2j*phase)
and do an abs() and 20log() from this value

I am getting incorrect values - instead of curve with max/min gain in -9 dB and steeper slope i have -18 dB.
Maybe i am using phase incorrectly?


How to sync audio and video signals in MATLAB? I have a video clip that has a frame rate of 25 fps (MJPEG compression) and an audio clip with a sample rate of 48kHz. Both of them play for the same duration - 27 seconds.

Is there any standard way to sync the two signals of different sampling rates if I wish to write them to file?

If this is not doable cleanly in MATLAB, I am happy with a some software tool on Linux or MacOS that can do this syncing.

E.g. the VideoFileWriter object can take both an audio sample and a video frame at a single step, however I have about 1917 audio samples corresponding to a single video frame.

Some details:

I have read in the audio signal using [signal, fs] = audioread(filename).
signal is a 1336321 x 1 array.

All the 697 frames in the video are saved as a cell array of RGB images.

Thanks to helpful comments from @MBaz, I managed to come up with a solution:

we can associate multiple audio samples with a single frame using the VideoFileWriter object. This fact and use-case is missing in the documentation.

First, some stats about the audio and video files. The stereo audio samples are in a 2xN array signal. The video frames are in a cell-array of RGB images, frames.

%% Write audio and video to file
%   Write both audio and video samples into a single video file.
%   Multiple audio samples are matched with one video frame.
%   MATLAB problem: Compression is not possible when audio is included.

% Display A/V stats
fprintf('\nVIDEO STATS:');
fprintf('\nNum. of frames: %d', length(frames));
fprintf('\nFrame rate: %.02f fps', annot.video.frame_rate);
fprintf('\n\nAUDIO STATS:');
fprintf('\nNum. of samples: %d', size(signal,1));
fprintf('\nSampling rate: %.02f Hz\n', infoAudio.SampleRate);

Here's the portion that matches multiple audio samples to a single video frame.
Note that there are 697 video frames and 1336321 audio samples. Some basic fiddling to ensure proper matching. The audio subsampling effect is imperceptible in the output AVI files that I tried.

numAudio = size(signal,1);
numRep = floor(numAudio/length(frames));
numDiff = numAudio - numRep*length(frames); % mismatch

if numDiff
    % if length(frames) does not evenly divide nAudioSamples, then
    % subsample audio to match numRep*length(frames)
    selector = round(linspace(1, numAudio, numRep*length(frames)));
    subSignal = signal(selector, :);
end
assert(numRep*length(frames) == size(subSignal,1));

Finally, we use the VideoFileWriter object available in the Computer Vision System Toolbox of MATLAB to write audio and video to file.

videoFWriter = vision.VideoFileWriter(fullfile(shotPath, 'avclip', ...
                                        [num2str(shotNum) '_av_clip.avi']), ...
                                      'AudioInputPort', true, ...
                                      'FrameRate',  annot.video.frame_rate);
for i = 1:length(frames)
   fprintf('Frame: %d/%d\n', i, length(frames));
   step(videoFWriter, frames{i}, subSignal(numRep*(i-1)+1:numRep*i,:));
end
release(videoFWriter);

U can use simulink
we can associate multiple audio samples with a single frame using the VideoFileWriter object. This fact and use-case is missing in the documentation.

First, some stats about the audio and video files. The stereo audio samples are in a 2xN array signal. The video frames are in a cell-array of RGB images, frames


Problems Using FFT to Compute Impedance in a Model Neuron I'm a neuroscientist currently investigating the resonance properties of a single neuron model that a colleague and I have constructed. The language we code in is Julia, which I hope is similar enough to Python for most people to understand. I've run in to a very particular problem, so bear with me please.

In order to investigate the resonance properties of the model neuron, I am applying an exponentially swept sinusoid current to it and analyzing its voltage response. In my field, such a current is called a ZAP current, though a more familiar term is probably Chirp signal. Just so there is no confusion, the current is generated from the following function \begin{gather*} ZAP(t) = I_{app}sin(2\pi f(t)t)\\f(t) = (0.001)\frac{f_{min}}{Lt}(e^{Lt}-1)\\ L = \frac{\ln{\frac{f_{max}}{f_{min}}}}{t_{max}} \end{gather*} where $t_{max}$ is the time duration of the simulation, $f_{min}$ and $f_{max}$ are the starting and ending frequencies of the sweep, and $I_{app}$ is the magnitude of the sinusoid. The voltage response of the model neuron to this input looks like this:

Clearly, there is a resonance in the model at whatever the input frequency is at around 50,000 ms. This means there should be a peak in the impedance of the model at this frequency. The impedance is calculated as the ratio of the discrete Fourier transform of the output voltage to the input current $$ Z = \frac{FFT(V(t))}{FFT(ZAP(t))} $$ I then plot the magnitude of the impedance against the natural logarithm of the frequency. I expect to get something like the one on the left, but instead I obtain the one on the right:

As I vary the parameter $E_K$, the resonance frequency changes, and I expect to see this reflected in the impedance-frequency plot, but every one of the impedance-frequency plots looks almost identical to the one shown above despite the voltage responses looking quite different.

I have spent hours going through my code line by line and experimenting with different modifications to figure out where I'm going wrong, but have not been able to produce the desired output. Here is the stripped down code:

## Time Parameters
dt      = 0.01         # ms
tmax    = 100000        # ms
tspan   = convert(Int64,tmax/dt)
t       = linspace(dt,tmax,tspan)

## Simulation and Analysis
# magnitude of sinusoid
Iapp = 1.0
# minimum frequency
f0 = 0.1
# maximum frequency
ff = 15.0
# create an array of the frequency sweep values for plotting
Zfreq = ZAPfreq(f0, ff)
# optimize Fast Fourier Transform for an array of size and type ZAP()
P = plan_fft(ZAP(Iapp, f0, ff))
# calculate FFT for input ZAP()
inFFT = P*ZAP(Iapp, f0, ff)
output    = simulate(ZAP(Iapp, f0, ff))
# calculate the FFT of the membrane voltage response
outFFT = P*output
# calculate the impedance of the cell at each frequency ZAPfreq[i] from Z = V/I
Z = outFFT ./ inFFT
# calculate the magnitude of the complex impedance
MagZ = abs(Z)
# center the DC response (f = 0 at center of frequency axis)
shiftZ = fftshift(MagZ)
# plot the impedance at each frequency
plot(Zfreq, shiftZ, grid = false,
  legend = false,
  linecolor= :black,
  xscale= :ln,
  xlabel = "log(\$f\$)",
  ylabel = "Impedance (Ω)",
  title = "Kispersky-Caplan Impedance Response with \$E_K\$ = $EK",
  top_margin=5*mm,
  bottom_margin=5*mm,
  left_margin=5*mm,
  right_margin=10*mm,)

Here, the $\texttt{simulate()}$ function is the one that actually carries out the simulation of the model. Any help in determining the problem is greatly appreciated.

I don't know how much background you have with FFT, so I apologize if I list a few things you already know.

In general, I think people have a tendency to throw an FFT at a problem and hope it works, but you may need to consider what the transform is really all about. First, when you take an FFT, it's implicit in the calculation that your signal is periodic. So, even though you're feeding it a finite number of samples with frequency content only between f0 and ff, the calculation makes the assumption that before time zero and after your samples, the waveform repeats forever. If your waveform's beginning and end don't match perfectly, there's going to be a lot of frequencies you wouldn't expect in your FFT. Basically, it's like there's a discontinuity. So, often, you'll window the signal to force it to zero at both ends so each periodic clone of your waveform meets the next smoothly. This obviously also has consequences for the frequency content of the FFT, so you need to do this carefully also.

Second, if your waveforms have any DC content (a constant offset, frequency = 0 Hz), it usually makes sense to remove it before your FFT unless DC is of interest to you. Often, low frequencies are huge and swamp the stuff you care about. Doing some high-pass filtering with a very low cutoff will help you see the interesting content of your plot more easily.

Third, when you view your plot on a linear frequency axis, you'll see that there are Fourier coefficients for both positive and negative frequencies, and the plot is symmetrical about 0 Hz. For qualitative assessments, you probably only need to look at the positive half.

Fourth, consider that you're taking the ratio of output over input, and anywhere your input FFT is near zero, that's going to blow up, and may make it hard to see what's going on.

So, I'd say, do some preprocessing of your signal to remove things you don't care about, and to prevent edge effects. Then post another plot and let's see how it looks. Since you're simulating the input and output here, there may not be a need, but in practice, I'll often use Welch's method to compute power spectral densities, and then compute an H1 or H2 estimation of the frequency response. Here, we can hopefully get a result with a simple analysis.

I have no idea about Neuron modelling so pardon me if my assumption is wrong but I can see that there appears to be a problem in in your electrical model.

Impedance is :  $ I(\Omega) = V(t)/I(t).$ But ratio of frequency volatge (V(t)) and current I(t) is not impedance. This quantity (ratio of frequencies), might have some significance in your model but that is not impedance. I would suggest trying FFT of the impedance function instead, which is :
                     $ Z = FFT(\ V(t)/ZAP(t)\ ) $

Also you seem to interested in input frequency at a particular time (where there is resonance), you can so this by STFT (Short Time Fourier Transform) or Wavelet transform.


Convolution of Two Impulse Signals I have encountered convolution of two different impulse signals.

x[n] = (1/2)^n . u[n-2] * u[n]
x[n] = u[n] * [n]

u[n] = discrete impulse signal
. = product operation
* = convolution operation

For the first one, I found this solution:

x[n] = 1/4 if n = 2
x[n] = 0 if n != 2

For the second one, I found impulse signal itself

Edit: Are my answers are true ? My professor told me that the answer for the first one is wrong, but he did not say the correct answer.

Your answer is right, assuming you posted the question right.

But you better use the standard notation as Dilip Sarwate already indicated; $u[n]$ is the unit-step and $\delta[n]$ is the unit impulse. Then

$$ 0.5^n \delta[n-2] \star \delta[n] = 0.5^2 \delta[n-2] = \begin{cases} { 0.25 ~~~, ~~~n= 2 \\ 0.00 ~~~,~~~n \neq 2 } \end{cases}
$$

you can get the answer for the second case, exactly in the same way.

This is basically an exercise to test the student's understanding of the concept that the unit impulse is effectively the unit in convolutions, that is:  $$\delta \star x = x$$ for all signals $x$. Perhaps a systems explaination might help. If an LTI system has impulse response $h$, then we know that the output of the system when $x$ is the input is $y = h \star x$. So, $\delta \star x$ can be thought of as the output of an LTI system with impulse response $\delta$ when the input to the LTI system is $x$.  What LTI system has output $\delta$ when its input is the unit impulse $\delta$?? It is just the canonical straight wire with (no) gain that audio enthusiasts dream about!  And so, $\delta \star x = x$ for all $x$.

With this, it it is easy to verify that the OP's answer to the first question is correct (but maybe his professor wanted to see the answer as $\left(\frac 12\right)^n \delta[n-2]$ or $\left(\frac 12\right)^2 \delta[n-2]$ instead of what the OP wrote) while the second answer $\delta[n]\star [n] = \delta[n]$ is incorrect, it should be $[n]$.


How to analyse periodic data obtained from a vibration experiment? I have data from a vibration experiment on a structure. This data was obtained by attaching an accelerometer to a FFT analyser, which records the data at a high sampling rate (typically in kHz range) and can also perform real-time FFT of the recorded data.

The image above shows the recorded signal below and the FFT above. The FFT shows all the dominant frequencies in the original signals.

What more patterns/ important data points can I discern from this data, that might have utility ?

For example, one being, the RMS value of the signal. What else ?

Also, the recorded data here as a mat file, which can be loaded in matlab, the data in a variable called Track1.

The physics of the problem: The data is collected from a shake table, a rigid platform, that oscillates back and forth, so the sharp edges in the original signal periodically, representing the point of return.

Remember that you are observing a system (the structure) driven by an input (shake table).  Assuming linearity, the FFT result was caused by the shake table and the structure.  You should be able to pick out the various components to identify mechanical resonances in the structure.  For instance, you should see a strong component at the frequency of the shake table.  Exactly how to do this in more detail is out of my domain, as I am an electrical, not mechanical, engineer.


Is interpolation of an audio signal to increase frequency resolution possible? I apologize if some of what I ask is not entirely correct, I'm new to this field, but extremely interested.

I have an Audio signal of sample rate 44.1 kHz that I want to segment into 30 frames, and get the DFT of each frame to find the magnitude of certain frequencies for that frame. However, this means that I have a frequency resolution of 30 Hz in each bin, which isn't narrow enough.

Is it possible to interpolate the data to attain more data samples? As far as I'm aware, doubling the number of points would give a 88.2 kHz sampling rate, but still give a frequency resolution of 30. Would it be possible to treat the interpolated data as still having a sample rate of 44.1 kHz?

I do not understand why you have 30Hz resolution so I will focus only to the principle of the question "does interpolation increase resolution?".

Short answer is no, no new data, no new information.

A longer answer needs the spectrum visualization below, with time domain, continuous frequency domain and discrete frequency domain from left to right.

The interpolation technique here is to preserve the information of spectrum. DFT works on discrete frequency domain which is the part from $-f_s$ to $f_s$ of the continuous frequency version.

First, look at the continuous frequency domain, if you upsample your signal correctly, it is equivalent to changing the sampling frequency. You wish to double the number of data, but it is just removing one-half the spectrum replicas of sampling process.

Now, look at the discrete frequency one. This version is normalized from $-f_s$ to $f_s$ of the continuous frequency counterpart. $f_s$ is doubled, the spectrum is then shrinked by a factor 1/2. If we call $0 < \alpha < 1$ the proportion of non-zero frequencies before upsampling, this proportion is $\alpha/2$ after upsampling. Before upsampling, DFT gives you $N$ bins for $\alpha$ then $N\alpha$ bins for the spectrum; after upsampling it is $2N$ for $\alpha/2$ then always $2N \times \alpha/2 = N\alpha$ bins for the same spectrum. No, your resolution does not change at all.

To have smoother resolution, the only way is to add more data. In your example, instead of dividing to 30 frames, divide your audio file to 15 frames.

Interpolation in the frequency domain does not increase peak separation resolution, but it can increase frequency estimation resolution if the signal to noise ratio is high enough.

You can interpolate the FFT results either by Sinc kernel interpolation (perhaps windowed), or by using a much longer zero-padded FFT, which both produce similar results.


SRRC filter -- filter span tradeoffs I am simulating (in Matlab) the transmit portion of a comm system.  I generate random data, then I want to modulate it to BPSK, pulse shape the data, and then potentially shift it up in frequency by a carrier.  My question is regarding the pulse shaping portion.

Atleast in matlab, to generate the root raised cosine filter using the rcosdesign function, you need to know the roll-off factor  ($\beta$) , as well as the samples-per-symbol ($sps$) and the filter span ($span$).  The samples per symbol will be dictated by the sampling frequency $f_s$ as well as the desired length of the generated signal with a certain datarate $R$.  For example, if I have a BPSK signal with a sampling frequency of 8192 Hz, and I want to generate 1 second of data at a datarate of 256 Hz, then

$$
sps = f_s/R
$$

Now we've established 2 of the inputs to creating the filter are set, the rolloff factor $\beta$ and the samples per symbol, $sps$.  Now, my understanding is that the filter span is kind of arbitrary.  I know that a longer filter span creates a longer group delay, but what are the other tradeoffs to having a longer filter span versus shorter?

Lower rolloff means less excess-bandwidth and therefore longer time span.
That is, you will need a longer filter to accurately reproduce the ideal RRC filter response (increases implementation complexity)
Longer filter means more accurate RRC filter and potentially lower ISI at the receiver side (depends on your target SNR)
Longer filter means less noise out of band transmitted (infinite length filter would have zero power outside symbol-rate*(1+rolloff)

In practice, due to the above constraints, RRC filters are usually not designed by just truncating the ideal RRC impulse response. Either you multiply the impulse response with a window function or you directly optimize the filter coefficients for error and rejection.

I'm assuming that you have got your roll-off from somewhere else Kiran, you say it is set but not how.  Not sure how the Matlab function defines filter span ( samples or symbols ) but you usualy get nicer spectra and eye diagrams if you choose a span where at least one of the end points is a zero.


Help in viewing the output of a neural network I used nftool in Matlab 2012 and trained a network. I gave the training inputs as x=[250:1] and targets as t=[250:1]. I used 10 hidden layers. I trained the network and got the results and exported the data to workspace. I finally got net as the variable in the workspace. I want to realize this network in hardware. So, I want the exact weights and bias values. I searched the net and typed net.IW and got 10 values, net.LW and got yet another 10 values and net.b and got 11 values. My network diagram is as shown below:

If IW is the Input Weight and LW is the Layer Weight, I realized the network diagram as shown below:

I got 10 values for IW, which would fit into 10 hidden layers, and 10 values for LW, which will fit into 10 hidden layers according to my realization. I got 11 values for bias, which will fit in to my network, since there are 11 (b) blocks in the network.
But, i am missing one IW and one LW value.

I want to know whether there is a mistake in my realization of the output of the network or I missed any of the values. Please help.

If I understand correctly your IW should be [1*10] as it takes 1-dim input to 10-dim hidden layer
than
LW should be [10*1] as it takes 10-dim hidden to 1-dim output,
b should be 10 for the hidden layer and 1 for the output layer

Why do you think your'e missing something?


How to Generate False Color Palette How to convert grayscale intensity input to something like the following image?

It seems that no response (0, black) is converted to dark blue, then it goes through other colors to red.

How are these colors derived? Is there any standard gray-to-RGB conversion function for these?

I have tried representing intensity as hue in CIE-LCH model and then converting to RGB through CIE-Lab. But the result weren't as nice as these MRI and PET scans.

Another approach may be a rainbow, i.e. going from high frequency (violet-blue) to lower frequencies (dark red). But is there a conversion formula for this?

What you're looking for is a color map. As its name implies, this process maps an input (single-component, or grayscale) value within a predetermined range to some other color based upon some mapping. The implementations I've seen typically use linear RGB interpolation between a number of control points specified throughout the allowed input range. The color map used in the image you showed is a pretty common one, starting with dark blue, then through cyan, green, yellow, orange, and ending with red. It is referred to as the jet colormap in MATLAB.

So, as an example, for an input intensity on the range $[0,1]$, you might have an RGB lookup table of:

0.0  -> (0, 0, 128)    (dark blue)
0.25 -> (0, 255, 0)    (green)
0.5  -> (255, 255, 0)  (yellow)
0.75 -> (255, 128, 0)  (orange)
1.0  -> (255, 0, 0)    (red)

Then, if you wanted to map the value $0.2$ to a color, then you would linearly interpolate between the first two datapoints as follows:

color = 1 / 5 * (0, 0, 128) / 0.25 + 4 / 5 * (0, 255, 0)
      = (0, 204, 26)

which would yield a greenish cyan, as it's close to the green control point.

Class to translate the example shown between [0, 1 ]

//Color Gradient [0, 1] => Result (R, G, B) R, G, B [0, 1]
void GradientCor::set_value (double y, int y_min, int y_max)
{
    /*
      y = mx+b
      x -> y
      -0 -> -80
       1 -> 120
      m = dy/dx= 120 - (-80) = 200; b = y + mx = -80; x = (y-b)/m
    */
    valor = (y-y_min)/(y_max - y_min);
}
    /*
    0.0  -> (0, 0, 128)    (dark blue)
    0.25 -> (0, 255, 0)    (green)
    0.5  -> (255, 255, 0)  (yellow)
    0.75 -> (255, 128, 0)  (orange)
    1.0  -> (255, 0, 0)    (red)
    x -> y para (dark blue)
    0 -> 128
    0.25 -> 0
    m = dy/dx= 1020; b = y + mx = 128; y = mx+b
    result = y/255
    */
double GradientCor::get_R()
{
    if (valor <= 0.25)
        return 0;
    else if (valor <= 0.5)
        return(double)(1020*valor - 255)/255;
    else return 1;
}

double GradientCor::get_G()
{
    if (valor <= 0.25)
        return(double)(1020*valor)/255;
    else if (valor <= 0.5)
        return 1;
    else
        return(double)(-510*valor + 510)/255;
}

double GradientCor::get_B()
{
    if (valor <= 0.25)
        return (double)(-512*valor + 128)/255;
    else return 0;
}


Can Morse Code be encoded using "silence" for dots and dashes? Morse code dots and dashes are often displayed as a sound, or light that is on for a short or long period of time.

What if the only way to communicate was by knocking.... (or similar encoding) that can only encode using extended (or short) durations of silence.

Is it acceptable to use morse code in this way (where short or long durations of silence are the dots and dashes)?

i used to be a ham radio person (like more than 40 years ago) and still know Morse code.  it's like learning to ride a bicycle; you don't forget.  between the elements ("dits" and "dahs") of a single character, the space of silence is supposed to be the same as the ON time of a single "dit".  call that a single unit of time.  so a string of dits would be at exactly twice the frequency as a stream of dahs.  but the space of silence between characters is supposed to be three units, the same length as a dah.

so, if you encode only the onsets of each element, whether it's a did or dah, you can't tell the difference between a lot of dahs in a row in the same character (like zero or "0") or the same number of dits in separate letters (like five adjacent "E" in sequence).  from just the onsets of the beeps, a single "0" (which is "dah-da-dah-dah-dah") sounds the same as "EEEEE" (which is "dit-dit-dit-dit-dit" at the same rate).  nor could you hear the difference from "EENE" ("dit-dit-dah-dit-dit" at the same relative onset times).

if you make the space between letters long enough that there can be no confusion between separated characters and the closer elements within a character, you can figure it out from just knocking.  but, given the element and character spacing set forth by the format of Morse Code, you cannot tell "0" from "EEEEE" from "EENE" with just knocks.

Of course you can. You can do anything you want to the signal as long as there was a way to transform it back.


information filter instead of kalman filter approach I read many sources about kalman filter, yet no about the other approach to filtering, where canonical parametrization instead of moments parametrization is used.
So I would like to learn on examples how to use this filter with information matrix and information vector (not mean and covariance from KF).
Can anyone help me with this, give some examples, webpages, or parts of implementation code to be a basis to learn?
Really appreciate!

Information filters are an application of a learning algorithm. One example of an information filter would be a recommendation system. On the other hand, a kalman filter is restricted to linear tracking. Its not learning per se.

If you are looking for machine learning/AI based approaches to SLAM, then you should probably look at Particle Filters.


comparing 2 audio files in c# I have 2 audio .wav files recorded by the the same person with the same word. Each audio file has a format of 16bit, 44.1kHz and is a mono channel. I would like to know how to compare these 2 audio files if it is by the same person. What I know is to process the files to produce a floating point result of values between -1 and 1 and apply fft on it and compare the frequency. What I would like to know is how to compare the audio files, prefably with the code which shows the comparison and if possible the steps before that(processing to become values btw 1 and -1).Thanks. Im relatively new in this area so I would appreciate If u guys would help me out. I'm currently using a metro application in vs. I would prefer not to use any library like naudio or bass.net but I would still accept any answer.
A step by step solution starting from utilizing the .wav file is much appreciated. I have been trying to search topics related to what I'm doing but what I have found are small parts of it.

I'm thinking the signals could be compared by

$D = \mathcal{F}^{-1}(\mathcal{F}(A) - \mathcal{F}(B))$

Where A and B are the two audio files to be compared.

I guess after that you could come up with some criteria for deciding that the voices are different. Perhaps when $|D|$ goes over some threshold.

You could also use Cross correlation

Using some other library to help with this would save time. But I haven't done this with audio files.

AForge Looks like it has FFT although I haven't used it.


FFT doesn't always have frequency in all frames I'm trying to do an FFT on live audio to see its spectrum.

So I have a sampling rate of 44100 Hz and from this take out 1024 bins on my FFT

So to test it I play a continuous tone at 15000 Hz and another one at 16000 Hz and look for it on my spectrum. I try to use higher ones to cut out noise from around. I have made the sound in a way there is no clipping

The thing is when I look at the individual frames of the FFT in each 2048 byte block of data I can see the peaks at these frequencies but some frames won't have the peaks, yet the others would.

How can this happen? How can I avoid it or what can I do to the way I'm reading the audio so that I can see the two tones in every frame?

You might be looking only at the cosine component (sometimes called the even or real component) of the FFT result, instead of the magnitude of each entire complex FFT result vector.  The real component will vary as your test signal goes in and out of phase with the cosine basis vectors (always even functions) of each frame, even if the complex bin magnitude stays constant (or varies only slightly due to windowing scalloping).

mag[i] = sqrt(real[i]*real[i] * imag[i]*imag[i]);

Not sure I understand your question exactly, but make sure to take the magnitude not just the real/imaginary components. Also, it will probably help to window your incoming data, and it sounds like you aren't doing that. For further detail and even some sample code, see here:

http://blog.bjornroche.com/2012/07/frequency-detection-using-fft-aka-pitch.html


Most appropriate spectral analysis method for a recording of a saxophone note? I'm currently doing an investigation where I am attempting to resonant frequencies of an Alto Saxophone at various different notes. I have taken audio recordings of the notes with a sample rate of $44100\textrm{ Hz}$, and exported them into CSV files. Now, I am attempting to analyze this data using the signal processing plugin for MATLAB, but I'm feeling a little bit overwhelmed with all the different options, as I'm not sure which spectral analysis method is most suitable for the my task.

In MATLAB, I'm seeing a number of spectral analysis methods including FFT, Burg Method, Welch Method, Multitaper Method, and more - but I'm not sure of their specific applications.

Again, my goal is to clearly identify the major spectral peaks (formants) which correspond to the resonant frequencies of the Saxophone bore. The data was taken with a condenser microphone at a sample rate of $44100\textrm{ Hz}$, with the Saxophone bore exactly 1 meter from the microphone. Ultimately, I'm looking to create a plot similar to the image I have attached.

Please excuse my ignorance on this subject, this is my first time using MATLAB and I am not well-versed in digital signal processing. If anyone could help me with this, I would be extremely grateful.

The plot OP shared looks similar to log magnitude spectrum(black line) plotted with its envelope(red line). so from what OP says, OP need to plot log magnitude spectrum of recorded audio/signal. One thing to remember for this analysis is that same note should be played for entire duration of the signal to be analysed. if not, one should use proper windowing methods to ensure that signal remains stationary(in same note) for the entire duration of window under analysis. I am assuming that a note corresponds to a few particular frequencies.

that is the way if you want to plot/analyse a graph similar to the one shared. are you also looking for other methods for determining the peaks corresponding to a particular note at a time?


Can we normalize data to unit variance without distorting the phase? Consider we have one dimension time samples. The noise in the samples made some portions of the time samples amplitude very high. Obviously the short variance for those time samples will have high variance. I want to do running variance (of some window size)  and normalize data by dividing by the variance. Is this procedure would destroy the phase of the phase of the signal? Is there a way to normalize data with unit variance without distorting the phase?

"Normalizing the variance": Since variance is a real number, normalizing it means multiplying the signal with a real number.

That doesn't change the phase, ever.


How does taking the absolute value of a complex signal reflect in the frequency domain? I have a frequency-domain representation $X(e^{i\omega})$ of the complex discrete one-dimensional signal $x[n]$: $X(e^{i\omega})=\mathcal{F}\{x[n]\}$. Is there a frequency-domain transformation of $X(e^{i\omega})$ into $\hat{X}(e^{i\omega})=\mathcal{F}(|x[n]|)$?

Obviously, $\hat{X}(e^{i\omega})=\mathcal{F}\left\{\left|\mathcal{F}^{-1}\{X(e^{i\omega})\}\right|\right\}$ will do the trick, but I am interested if there is a dual to the absolute value function $|x|$ in the frequency domain.

Absolute value is not a linear operation, so no, there is not a consistent dual in the frequency plane.

Having said that, though, it is somewhat similar to calculating the signal power, which involves multiplying $x[n]$ with it's complex conjugate.  Multiplication in the time domain is equivalent to convolution in the frequency domain.  The result in the time domain is real, so the result in the frequency domain is even (i.e. symmetric around 0 Hz).

Absolute value will probably produce a similar effect to calculating the power since it is the square root of the power.  It is also real so it will be even in the frequency plane.

Well the Fourier transform of the absolute value function exists:

$-\frac{\sqrt{\frac{2}{\pi
   }}}{w^2}$

Source:
FourierTransform[Abs[x], x, w]


High pass filter of 0.01 Hz for EEG slow drifts I am trying to make a high pass filter for electroencephalographic (EEG) datasets, in order to get rid of very slow drifts. However frequencies around 0.3 Hz are very important for the research these data serve.

I use Matlab and I tried various things. It seems that a very low dB attenuation, of 3 dB is fine with FIR filters. For example I used that:

h=fdesign.highpass('Fst,Fp,Ast,Ap',0.005,0.01,3,1,250); % i needed to cut
d=design(h,'butter');
fvtool(d)

However this filter makes a dc drift at the start of the data and I would not wish to loose these first datapoints.

I have read through your pages that a clever idea is to design a low pass filter that later I would subtract from the real data. I did it with an equiripple filter and it worked but left a baseline lift. I do not want to run now average removal as this is in a next step in the analysis protocol.

Any suggestions?

You're asking for a time machine: in order to avoid the DC bump at $t=0$, you need to know the state of the filter as if it had been running before the recording ever began. That's not going to happen no matter what kind of filter you use.

Here's a trick which may come in handy for you. Let's say you have an $N$-tap FIR filter (or an IIR filter whose impulse response decays adequately to 0 after $N$ samples). Take the first $N$ samples of your signal, reverse them, and prepend that to the start of your signal. In effect, we're defining a new signal $g(t)=g(-t)=f(t)$ for the purposes of filtering. When you filter that, the filter state will have already been "primed" with the initial signal: delete the first $N$ samples of the output, and the DC bump should have hopefully gone away.

Needless to say, there are numerous gotchas with that scheme — for starters, it's harder to do for realtime processing, and if the signal derivative is high at $t=0$, you are still going to have trust issues with the first few samples. But it should get rid of the bump.

You don't care about absolute DC value, correct?

Why not just add a static offset to your data so that the first data point is at zero?

Sure, you would have to add the offset to every datapoint, but it would completely avoid the big step-response issues you are having with the filter.

Effectively the filter starts up initialized to zero. Therefore, when your data begins, the filter sees a big stair-step from zero to whatever DC level you have.

Just add offset to remove the stair-step.


Audio sounds oddly distorted https://clips.twitch.tv/TangentialObliqueLobsterSaltBae
https://clips.twitch.tv/CharmingBashfulAlfalfaRickroll
https://clips.twitch.tv/NaiveKindGorillaGivePLZ

(Disclaimer: These were not recorded on or uploaded from my computer.)

In the 1st video linked above, the water sloshing sounds heavily distorted and the same type of "noise" heard throughout that segment of the clip can be heard periodically throughout this person's videos, including the other ones I linked.

I've been trying to help this person hear the problem with his audio, but he's not hearing it and I'm not entirely sure what the cause of it could be or how one could even describe it.

EDIT: He switched to different broadcasting software and finally hears how bad it sounds now. Now we just need to know the cause.

It's mono, not stereo like one would expect, and it sounds like exacerbated compression artifacts to me, like in the difference signal between the stereo channels of a lossy compressed audio, but I'm not sure if that explains everything.

Given that he uses an external mixer for recording, there may be a simple explanation: He has connected a stereo cable to a balanced input of the mixer. The balanced input will typically take the left channel as "hot" and right channel as "cold", and will effectively see the mono difference signal: hot minus cold. He should use a stereo to 2x mono splitter cable, or use a stereo input in the mixer if available.

We determined that it was his broadcasting software that was causing the issue. He switched to a new broadcasting software and the audio sounds clean again. Thank you all for helping!


Transforming FFT to Detect Baud Rate I've seen in practice taking an FFT of a signal, then transforming it by the 2nd, 4th, etc. power (depending on signal type) to determine the baud rate. I'm trying to implement this in my code, but am lost as to what a transform of the 2nd power is. It appears to be some sort of non linear transform and not simply squaring of the FFT values.

As noted in the following link, baud rate is estimated from "the spectrum of the signal after nonlinear processing."  In their case, as noted in the first three sentences of the section entitled "Proposed Spectral Approach for Bit Time Estimation" - p. 916), the nonlinear processing is the first derivative followed by an absolute value operation:

http://www.iasj.net/iasj?func=fulltext&aId=51991

Others have proposed/analyzed estimation using different nonlinear operations, such as squaring of the signal or using higher derivatives.  Obtaining the references for the above paper and doing a thorough search would be a good place to start to learn more.

I'll bet that the idea was roughly as follows:

Baud rate detection is done via the autocorrelation function of the input
A correlation is expensive to perform, much more than multiplication.
A correlation in the time domain is equal to a multiplication in the Fourier domain
An auto-correlation in the time domain is equal to squaring in the Fourier domain

But remember that the FT is complex-valued. Ignoring the phase from the FFT will cause unpredictable results after the IFFT. That the phase into account, and the result of your IFFT should be real again (minus rounding errors)