Creating derivative images, based on original, that can only be visually understood when stacked I've been thinking about this for a while. What I have been trying to accomplish is a way of taking an image (png specifically, or any other format that allows transparency) and generating N derivative images. Each derivative image would, theoretically, be as difficult as possible to visually connect to the original. If I processed an image of a mountain into three images, none of the three images should suggest to you that it is a mountain. However, if all three images were stacked on top of each other, what you should see, with as minimal corruption or loss of detail as possible, is the original image.

Different approaches I've considered are splitting an image by color, and stacking the different colors. This works okay when N is a large number. It works poorly when N is small. Quantization helps, but can be tremendously lossy at small Ns, where it would help most.

Replacing (N-1)/N pixels with noise or transparency in each layer works well at higher levels of N, but a bit poorly at N=2 or 3 because the top layer reveals the subject of the image.

It seems that dealing with the top-most image is where things become really difficult. I'm not really sure how to approach this.

...a way of taking an image (png specifically, or any other format that allows transparency) and generating N derivative images.

You can split the image into different spatial frequencies and recombine them to form the original. This is exactly the way that the Discrete Fourier Transform works.

Taking as a starting point the following image (found by Google Images and finally obtained from this article):

If you progressively filter it into different spatial frequencies, you will get intermediate images that look something like this:

Low Spatial Frequencies:

Higher Spatial Frequencies:

(...but not including the low, i.e. a bandpass filter)

Even Higher Spatial Frequencies:

(...another bandpass filter but shifted higher in the spatial frequency spectrum).

All the way to 11! Spatial Frequencies:

(...similar as above, this can be a bandpass filter or a high pass to include all remaining spatial frequencies).

Now, in your favourite image manipulation software (I am using GIMP with GMIC here), turn on transparency and put the different spatial frequency bands in different layers so that you can turn them on and off. You get something that (progressively) looks like this:

All that you need now, is some transparency "paper" and a laser printer. And also to note that if you mount the transparencies one behind the other, then you might have to adjust the scale of each image depending on the distance they will be mounted at. You can, for example, mount one on your door window and the other on the window behind the door so that when you look at the images from a specific point, they (both) make sense but the one that is further away will have to have its size slightly enlarged.

You will of course have to experiment a little bit with different spatial frequencies and different spatial frequency bands to adjust the images to the result you are trying to achieve. I only used a few bands here.

Hope this helps.

EDIT:

In the case of limited dynamic range (e.g. pixels can only be on or off), the technique remains the same but the gradient will have to be remapped through the use of halftoning.


What is the effective alpha of 2 cascaded EWMA IIR filters? Considering simple 1st order exponentially weighed moving average (EWMA) filters of form

$$y[n] = \alpha x[n] + (1-\alpha) y[n-1]$$

what is the effective alpha of 2 cascaded EWMA IIR filters?

Maybe I'm asking the wrong question, but not having a lot of luck Googling for the answer

It is important to define in what sense two cascaded EWMA filters should be "equivalent" to a single EWMA filter. One possibility would be to define an "effective memory" or "effective averaging window size" and say that the two filters are equivalent if they have that same effective memory.

The second question to be answered is how to define the effective memory length of such a filter. Again we have to make a reasonable choice, and one possibility is to look at the filter's impulse response $h[n]$ and choose the length $N$ such that the sum over the first $N$ coefficients equals a certain percentage of the sum over all coefficients:

$$\frac{\sum_{n=0}^{N-1}h[n]}{\sum_{n=0}^{\infty}h[n]}=r,\quad 0\ll r<1\tag{1}$$

That is, we consider the coefficients with indices greater than $N$ negligible. The impulse response coefficients $h[n]$ are the weights to compute the output signal from the input signal:

$$y[n]=\sum_{k=0}^{\infty}h[k]x[n-k]\tag{2}$$

For the filter given in your question it can be shown that the impulse response is given by

$$h[n]=\alpha(1-\alpha)^n,\quad n\ge 0\tag{3}$$

Furthermore,

$$\sum_{n=0}^{\infty}h[n]=1\tag{4}$$

and

$$\sum_{n=0}^{N-1}h[n]=1-(1-\alpha)^N\tag{5}$$

With $(4)$ and $(5)$, Eq. $(1)$ reduces to

$$1-(1-\alpha)^N=r\tag{6}$$

which, for a given percentage $r$, can be solved for either $\alpha$ or $N$:

$$\begin{align}\alpha&=1-(1-r)^{1/N}\\
N&=\frac{\log(1-r)}{\log(1-\alpha)}\end{align}\tag{7}$$

So, e.g., with $r=0.9$ and $\alpha=0.1$, the effective averaging length is $N\approx 22$. Alternatively, for a desired length $N=100$, $\alpha$ must be chosen as $\alpha=0.023$.

The same calculations can be done for the cascade of two such filters, the only problem being that it is not possible to give analytical expressions analogous to Eq. $(7)$. However, the finite sum given in $(5)$ can be easily computed using the filter itself. If you cascade two EWMA filters and use a constant input $x[n]=1$, $n\ge 0$, the output will be the sum $(5)$. The first output sample is the sum for $N=1$, the next one for $N=2$, etc. So by this computation one can figure out the value of $N$ for a given percentage $r$. This value is given by the number of output samples necessary for the output to exceed the given value $r$.

It's time for an example. Let's again choose $r=0.9$, i.e. the effective memory size is define by the number of coefficients which add up to $90\%$ of the total sum of all coefficients. Let the two $\alpha$-values of the cascaded filters be $\alpha_1=0.1$ and $\alpha_2=0.2$. To find out its effective averaging length (at least according to our definition), we feed the filter with a constant input $x[n]=1$ until the output exceeds the value of $r=0.9$. The output is shown in the figure below. The output values are blue, and $r$ is shown by the red line.

The effective averaging length $N$ turns out to be $N=27$. In order to find out which first order averaging filter has the same effective length, we need to compute $\alpha$ from Eq. $(7)$ for $N=27$ and $r=0.9$, which gives $\alpha=0.082$. Hence, in the above defined sense, we can say that two cascaded EWMA filters with $\alpha_1=0.1$ and $\alpha_2=0.2$ are equivalent to a single EWMA filter with $\alpha=0.082$.


Inter-Symbol Interference (ISI) channels Is there any open and free data-set for digital ISI channels? I want to have a variety of ISI channels' taps, such as $h[k]$ here. Especially looking for BPSK channels.

Things like the good ole' GSM standard specify power delay profiles and a phase variation model. So, there you go, reference model delivered in the shape of a standard.

Other than that, if you're just looking for any ISI-inducing (i.e. frequency-selective) channel, why not use your favourite DSP framework to come up with a random FIR filter?

Also, a lot of frameworks do come with their own channel simulations; GNU Radio for example brings multiple channel models of varying realism in its gr-channel module.

I often find that students overestimate the complexity of actually measuring a channel impulse response to use it in their own channel simulations. Just send a PN, receive it at the same time, and use a correlation to find the impulse response; if your receiver and transmitter exhibit a frequency offset, you might need to try a couple of frequency shifts and use the one that gives the highest correlation coefficients.


cyclic prefix in ofdm In ofdm systems after applying IDFT we take the last N/x samples and place them before these N samples. This is called cyclic prefix(CP).So why to use a cyclic prefix?

Why is the CP known to combat the  inter symbol interference.Why can't we use some junk data or zeros over there? Can't this junk data combat the ISI ?.

Regards,
Phani Tej

The cyclic prefix is to make sure that any multi-path interference (or other process similar to a linear time-invariant filter in the transmission channel) acts as a circular convolution on the FFT data frame, thus not affecting the orthogonality of the data channels within the FFT bins.  It also allow some slop in the receiver's symbol clock.

Thus the cyclic prefix has to be a circular repetition of the IFFT data, either before or after the IFFT block, and long enough to be as long or longer than the maximum expected multi-path delay time differential in the channel, plus the maximum of any timing error between the received data block and the receiver symbol clock.  The N/4 is just a default estimate of this time for some common channels, but longer may be required.

Anything other that a circular repetition can cause multi-path delays to act as a non-linear spectral spreading function (relative to the FFT frame) which will affect channel orthogonality, and thus S/N related reliability.

Repetition (circular or otherwise) also helps to statistically even out the PAP (peak to average power) ratio over the extended frame.  And prefix repetition may help push any transmitter anomalies due to sudden changes in this ratio forward of the FFT frame.


response of algorithm in non-equidistant time We are investigating ways to test a control algorithm.  The algorithm has a non-equidistant track of input data (i.e. not every sample is valid, and we know it), and should output a series of correction factors.

Under a certain frequency, it should correct for a deteriorating light source.

Is there a way to black-box test the algorithm's transfer function?  I'm having trouble especially with the fact that it's non-equidistant.

The goal of the tests would be to have a clear idea on the algorithm's behavior in terms of

delay/phase response
resonances
frequency response

So I was actually dreaming of a kind of a 'swept sine' test, or an impulse response, or something alike, but I'm a bit lost in what is feasible with non-equidistant signals.

I would try to stick with a conventional discrete time control approach.  Trying to deal with missing samples in the control algorithm complicates things and without dealing with them the control output will not be very stable.  Can something like this work in your application?

Here the assumption is that you can detect missing samples and use an interpolation algorithm to fill in the gaps before applying the signal to your controller/compensator.  The "Control Signal" input is a reference that can be used to set the normal operating intensity of the Light Source.

I would characterize the transfer function with using normal control signals (impulse, step, swept sine, all with uniform sampling).  Then I would model the real signal (from the lamp transducer) with some impulse variations that match the characteristics of your missing samples.  By impulse variations I mean adding small error quatities to signal to represent error from your interpolator.  Test how well the transfer function regulates in the presense of this "impulse noise".


Measure 3-D length using smart-phone's accelerometer I'm wondering if we can measure 3-D length only by our smart-phone's accelerometer. And we all know these low cost IMUs are not accurate.

You can model accelerometer's error this way:
$$ a = f*a' + g + b + \eta $$

where $a$ is the actual acceleration, $f$ is a 3x3 matrix modeling scaling, misalignments, cross-axis and ... errors. $g$ is gravity, $b$ is bias and $\eta$ is noise.

I'v found a good method to calibrate(compute $f$ and $b$) smart-phone's accelerometer:
Estimate smartphone accelerometer bias

Now I have to:

1) remove $g$

2) remove $\eta$

How can I do that accurately? As I mentioned above, this problem requires good accuracy.

If you believe accelerometers are not accurate enough for this problem, and other sensors should be used (Sensor Fusion), please share your idea.

Your model may be misleading, as bias and noise (including non-linearity and quantization noise) get double integrated, thus error increases quadratically and without bound over time.


Question about image equalization Let $f$ be an $N \times M$-image and suppose that $g$ is its  equalization.
Question  - Why the equalization does not change the objects, why a car on $f$ does not convert into a bus on $g$?. It is not obious that after the pixel intensity $r_k \mapsto s_k$ recalculation for $g$
$$
s_k =\frac{L-1}{NM}\sum_{j=0}^k n_k, k=0,1,\ldots,L-1,
$$
we obtain the same scene as was for $f$! Here  $n_k$ denotes the number of
pixels that have intensity $r_k.$

It probably has the potential to do so, but the probability is marginal. Result of equalization depends on the initial pixel intensity distribution.

Image histogram equalization can have different implementations.
In one possible implementation, you would re-map a whole block of a given intensity level into a new calculated intensity. In such a case, objects will not change shape; boundaries will remain the same.

The aim of equalization is not necessarily the creation of a new image with strictly flat histogram will all bins equally occupied. That won't create a good looking image in any sense. Instead, the aim is to expand the pixel intensities to fully utilize the available dynamic range in the hope of getting a higher constrast with more visible dark regions and unsaturated highlights.


Median filter algorithm in opencv I was wondering about the implementation of the  Median Filter in opencv because I am suffering form its slowness. So, I intend to implement another one on my own but first I want to make sure about the used algorithm in opencv.
Many thanks

I don't know about opencv's implementation, but a median filter can indeed be very slow because for every windows you have to sort your data in order to find the median. However there is a recursive implementation of the median filter, which is much more quick. It is very simple to build.

The basic idea is to recursively build the data's histogram that will allow you to very quickly get the median without sorting :
For the first window, fully build the histogram, then, shift your windows by one step (or more, it depends on you), remove the last value from the previous window and add the new one, you updated your histogram with only 2 operations.
Then you have to find the median within the histogram, this is recursively feasible but I don't remember the exact implementation (it is a bit more tricky), however simply checking the histogram's cumulative sum, while brutal, will give you the median quickly (worst case scenario is N summation, with N your window's size).


How is GPS used to synchronize oscillators in a distributed system? I heard that it is important to have all oscillator clocks at different places be synchronized for example in a geographically distributed system such a telecommunication network. Of course to some error margin.

They fix GPS receivers in the base stations to synchronize the oscillators. I know GPS transmits time, as in our watches.

Could someone please explain what does this time of the day has to do with synchronizing oscillators?

Thanks a lot.

Part of the GPS signal is a 1 pulse per second clock that can be used in distributed systems for timing synchronization.

This patent has some diagrams that might help, e.g.

Regarding how accurate it is, I'm not sure. This paper has the following graph... which seems pretty good to me.


Why is the median of a PDF not always equal to the mean? When studying probability density functions(PDFs), the focus is often placed on the "mean" and not the "median". The mean is often described as the center of mass of the PDF. I've always taken this as true.

However, the definition of the median seems to correspond to what I would consider the center of mass. That definition being:

$\int_{-\infty}^{median}pdf(x)dx = \frac{1}{2}$

Looking at this, we're essentially integrating the PDF until we reach the half-way point. To me, this seems identical to what we would consider the center of mass.

So why does integrating a PDF until we reach halfway (the median) not always equal the mean if the mean is the center of mass?

Note that I understand the difference between the mean and median from a general statistics perspective, but I'm curious from a PDF perspective.

I'll also format the question in the following manner:

$\int_{-\infty}^{\infty}pdf(x)dx = mass$

$\int_{-\infty}^{\infty}x*pdf(x)dx = center\:of\:mass$

If we define the following equation:

$\int_{-\infty}^{T}pdf(x)dx = \frac{mass}{2}$

Why is T not always the center of mass? It appears T is the center of mass if the function is symmetric, but why doesn't it hold true if the function is non-symmetric?

If you consider the pdf as defining the distribution of mass in a thin rod of infinite length and variable density (and so the "weight per inch" of the rod is different depending on which part of the rod we are talking about), then the mean is the point about which the rod is perfectly balanced (that is, if we set up a fulcrum and place the rod on it, it will be horizontal and balanced perfectly if the fulcrum is at the mean, and will tilt downwards on one side or the other if the fulcrum is anywhere else.

The median, on the other hand, is the point for which we can say that half the mass is to the left of the median and half the mass is to the right of the median. If that is what floats your boat as the one true and the only meaning of "center of mass" and you feel that everyone else who believes otherwise is in a state of mortal sin, so be it. Just be aware that when you refer to the median as the center of mass, you are likely to be misunderstood by those lesser mortals who have not attained your state of enlightenment and who still believe that center of mass is another name for what they call the mean of the pdf or random variable.

With respect to the comments on the main question as to an asymmetric pdf means that the mean and median will be different, this answer over on stats.SE (where this question more properly belongs) gives an example of an asymmetric pdf for which the mean, median and mode all have the same value.


How to calculate auto-correlation of a bpsk modulated signal or how to calculate expectation value of complex exponential function How to calculate auto-correlation of a bpsk modulated signal, or how to calculate expectation value of complex exponential function manually not by using matlab or any other software?

For example, if

$$x(t)= a\exp(i2\pi ft)+n(t)$$

how would you solve

$$\mathrm{E}[x(t)]=?$$

and

$$R(t,t_1)= \mathrm{E}[x(t+t_1/2)x'(t-t_1/2)]=?$$

where $a$ is either $+1$ or $-1$ mapped from the binary representation of the ASCII value of any text message.

This question is essentially unanswerable because too many details are left out.
Here is an answer to the following simplified version of the question

  What is the autocorrelation function of the random process
  $$X(t) = A\exp(j2\pi f_ct), -\infty < t < \infty$$
  where $A$ is a random variable equally likely to have values
  $+1$ and $-1$, and the autocorrelation function is defined as
  $$R_X(t,\tau) = E[X(t+\tau/2)X^*(t-\tau/2)], -\infty < t,\tau < \infty.$$

This is easy since the only random variable involved is $A$, and
it is a discrete random variable.  Thus, we have
$$\begin{align}
E[X(t)]&= E[A]\cdot \exp(j2\pi f_ct)\\
&= 0 ~~\scriptstyle{\text{since }} E[A]=0.\\
R_X(t,\tau) &= E[X(t+\tau/2)X^*(t-\tau/2)],\\
&= E[A^2]\cdot \exp\left(j2\pi f_c(t+\tau/2)-j2\pi f_c(t-\tau/2)\right)\\
&= \exp(j2\pi f_c\tau)~~\scriptstyle{\text{since }} E[A^2]=1.
\end{align}$$

What would be needed to make this question worth answering? First,
one usually has a baseband pulse shape, usually time-limited,
that is nowhere mentioned by the OP. The autocorrelation function
of a BPSK signal is usually a scaled version of the
autocorrelation of the baseband pulse shape (and so the power
spectral density is determined by the shape of the baseband pulse).
To get to this result, one needs to consider not just
one bit but rather the stream of bits that modulate the BPSK signal,
and assume that these bits are i.i.d. random variables (such as
$A$ above). Second, something needs to be said about the noise
$n(t)$ that the OP casually includes in her question.


How to determine best fit from two? I'm studying some variances ($\sigma^2$), that in my case, they must depend by velocity squared $v^2$.

So, in my experimental proofs, I have plotted a graph of $\sigma^2/v^2$, and I was hoping to find a horizontal line.

But, my data seem to be constants in $\sigma^2/v$, and now I find myself in times of trouble: I have to choose which one is the best fit.

In each cases ($\sigma^2/v^2$ or $\sigma^2/v$), I obtain a nearly straight line.

So, I look the equation of best fit lines as $y = mx +q$, where $m$ is his angular coefficient.

The angular coefficient $m$ is less, in the case of $\sigma^2/v^2$ (by about an order of magnitude), but in the case of $\sigma^2/v$, the range of the best fit line (from its maximum and minimum in the graph) is relatively more lower.

At this point, to choose the best fit, the best solution is:

Choose the case where $m$ is lower.
Or I have to choose the dynamic case (maximum value minus minimum value), where the best fit line is more flat in the graph (not by absolute value)?

Here are an example: in the first graph $m$ is greater (see the equation in the plot), but the line varies of 0.01.

In the second case $m$ is lower, but the line varies of 0.15.

What is the best fit? The first or the second?

EDIT

For Matt L. - Your approximation error $e$ is very similar to the $R^2$ parameter, that I define so (now I use your notation):

$$R^2 = 1 - \frac{\sum_i \left[d_i-f(v_i)\right]^2}{\sum_i \left[d_i-\overline{d_i}\right]^2}$$,

where $\overline{d_i}$ is the average value of $d_i$.

Maybe you would like say that I must simply see where I have the minimum error $e$ or the best fit parameter, e.g. $R^2$.

It looks like a straight line is not a great approximation in either case. One thing which is unclear to me is that the y-axis in the top plot is $\sigma^2/v$ whereas in the second plot it is just $\sigma^2$ (as far as I can see), and not $\sigma^2/v^2$ as I expected.

What you need to figure out is what you should actually expect. When you say that $\sigma^2$ should depend on $v^2$ then this doesn't say much. From what you do it looks like that you expect those two quantities to be proportional:

$$\sigma^2=\alpha v^2$$

because then you would ideally have a constant ratio $\sigma^2/v^2$. Is this the case?

Note that if you model $\sigma^2/v$ as a linear function of $v$, then you have

$$\sigma^2/v=mv+q$$

which is equivalent to

$$\sigma^2=mv^2+qv$$

So you actually model $\sigma^2$ by a quadratic function of the variable $v$.

To determine the best fit, it seems natural to compare the normalized approximation errors of the different fits. Given that $d_i$ are your measured data points, and $v_i$ are the discrete velocities corresponding to the data $d_i$, the approximation error can be defined as

$$\epsilon=\frac{\sum_i \left[d_i-f(v_i)\right]^2}{\sum_id_i^2}$$

where $f(v)$ is the approximation function, which could be linear as you've suggested in your question, i.e. $f(v)=mv+q$.


What are some ways to detect the fundamental frequency in a signal using without using FFT? I am trying to estimate the fundamental frequency in a signal to corroborate the result which I have got by FPGA implementation of an FFT algorithm. Hence, could someone give me some idea how to do so?

Count the time interval between zero-crossings (if the waveform/SNR allows...).
Autocorrelation, average magnitude difference function and their variants (YIN...).
Estimating the parameters of a damped sinusoids model (MUSIC/ESPRIT).
Several instances of Goertzel algorithm, or more generally other filterbanks (band-pass, comb) - in case there is only a small number of candidate frequencies that you want to discriminate.


converting block diagrams with delays to equations given this block diagram:

can someone please clarify how to interpret the second delay (the one prior to $Y$)? to convert this block diagram into a simple equation, I interpret the diagram as:

bottom path: take input at $n$, multiply by $-1$ and delay by $1$, gives: $-1x[n-1]$
top path: take input at $n$ and sum it by result of bottom path, gives: $x[n] + -1x[n-1]$
take sum of bottom and top path and delay it (?) - this is the part I'm confused about. it doesn't make sense to write: $(x[n] + -1x[n-1])[n-1]$ which is what the diagram looks like to me.

what is the right way to read this diagram?

If you are confused about the interpretation, why not divide the diagram into 2 parts. The signal to the left of the rightmost delay element is assumed to be z[n].

Thus, now $$y[n] = z[n-1]$$

$$z[n] = x[n] + (-1)\cdot x[n-1]$$

Thus

$$y[n] = x[n-1]  -x[n-2]$$


DSP finite impulse response (MATLAB) I designed the filter below on paper:

FIR
low pass
type 1
length 3
alpha = 6
beta = -12
allow freq. <= 0.1rad/sample
blocked freq. >= 0.4rad/sample

To confirm that my calculations are correct I was planning to design the filter on MATLAB. Is there a function to design this filter, apply an input frequency and check the output result?

Like Matt L. pointed out, alpha and beta is unclear as is 'type 1'.

You can use the command fdatool to interactively design filters in MATLAB, although you may need the DSP system toolbox installed.

From there, you can apply the filter using the filter() function.

As an example, if you design the filter using FDAtool, and the filter is an FIR, you can use filter(Hd.numerator,1,X,Fs) if Hd is your filter object, X is the input sequence, and Fs is your sampling rate. This is not the only way to do this, but just the way I would go about it.

X could be generated as a sin(2*pinf*t) where t = 1/fs and n is an index in a for loop or shortened to 1:ns for the number ns you wish to have your input sequence have as a length.


Random function covariance I was studying signal processing and I was frequently asked to verify if a certain covariance is possible for a given random function. I tried to check by verifying the property to prove it:
${γ}_{xx}(τ) = {γ}_{xx}(t+τ,t)$
but didn't quite got anything.
For exemple, if we suppose $τ_0$ a given constant, how would I prove that I can have the following covariance (or not)?
${γ}_{xx}(τ) = exp(-(τ- τ_0)^2)$

You have two different meanings for ${γ}_{xx}$ when you write stuff such as
$${γ}_{xx}(τ) = {γ}_{xx}(t+τ,t)$$
because on the LHS, the function ${γ}_{xx}$ has one argument while on the RHS it has two arguments. Be that as it may, a very important property of what you call the covariance but what most everyone else on dsp.SE would call the autocovariance function is that
$$\text{For all values of }\tau, |{γ}_{xx}(τ)| \leq {γ}_{xx}(0).\tag{1}$$
It is easy to verify that $\exp\left(-(\tau-t_0)^2\right)$ has value $\exp(-t_0^2)<1$ at $\tau=0$, and larger value $1$ at $\tau=t_0$, in contradiction of $(1)$, and so $\exp\left(-(\tau-t_0)^2\right)$ cannot possibly be a valid autocovariance (except, of course, for the special case of $t_0=0$).


Spectrogram of an executable binary One of the approaches I have been searching so far is the application of utilizing Neural Network to classify Spectrogram images, which is heavily used in speech recognition.

I am conducting an experiment to identify the similarities of binary executable (raw 0s and 1s), and planning to deploy the same concept of speech recognition on binary executable. I have been preparing the executable with the following criteria:

1- The sampling rate is 8 KHz.

2- The bit Depth is 8-Bits. This means that each sample's amplitude will be represented using 1 Byte.

I am being stuck in generating the Spectrogram. Now my background is networks and programming and I have been doing a lot of research to generate it. The challenge is the identification of the frequencies of the samples I have in order to have the 3-dimensions of spectrogram (frequency, time and amplitude). I did a research, where FFT can be used to generate the frequency domain from the time domain, but not sure if it is feasible for binary executable.

Python supports .wave files spectrogram plotting, but the files I have are executable and not music or sound.

1- Can I identify the used frequencies in my binary by having the sampling rate and amplitude?

2- Can I generate a spectrogram of an executable binary, even if it is NOT originally speech or music?

3- Is there any tool/code that I can utilize for the spectrogram generation?

This serves as potential starting points only, on how to transpose classical signal processing tools to specific kinds of binary data.

First, remember that standard spectrograms, fed as inputs to classifiers, can be seen as means to extract attributes from non-stationary data. There is apparently little need for inversion, so we are interested in analysis only (no synthesis).

Here, we have sequences of 8-bit words $w_0w_1w_2\ldots w_{N-1}$. When dealing with binary patterns, the frequency domain might not be the most appropriate. In that case, as with communication systems, the related notion of sequency could be useful. It relates to blocky oscillating patterns of $\pm 1$ (or $0/1$):

The number of zero-crossing increases in a manner similar to augmenting frequencies in Fourier transforms. I would thus prefer to use those vectors in sliding windows to analyze your binaries. You can find literature in such functions as Walsh, Paley or Hadamard (or WaLeyMard) sequences, depending on their ordering. The proper choice of the transform length is not evident.  This seems to exist already, under the name: "Walsh spectrogram": On Generating Walsh Spectrograms, 1976.

I would think it might depend a bit (pun intended) on the compiling or storage system. There is a possibility to look at dyadic size factors with Haar wavelets or Haar wavelet packets.

Note that there exist transformations for data taking values in finite fields, from The Fast Fourier Transform in a Finite Field, 1971, J. M. Pollard.

Other references:

Harmut, H. F., 1970: Transmission of information by orthogonal functions
Agaian, S. S., 1984: Hadamard Matrices and Their Applications


Processing delay window-function In some DSP slides, I've noted that when we try to choose a window for truncating a signal, sometimes a window is chosen that minimizes the processing delay (in my example there were a rectangular, Hamming and Bartlet)...but what is exactly this processing delay?

The minimum process delay I would have to have with Rectangular window?

Sorry, but why is Bartlet? I thought that the minimum phase was proportional at the width of the main lobe (I've resolved thanks, I did not remember that energy end amplitude is related.)

One interpretation of processing delay is that the causal window, $w[n]$ has more energy close to $n=0$ than any other window, $w_o[n]$ with the same frequency response (see section 5.8.4).

So that
$$
\sum_{n=0}^M |w_o[n]|^2 \le \sum_{n=0}^M |w[n]|^2
$$
for all integer $M \ge 0$ .

So if I plot these sums for the three windows, you can see that the Hamming window (barely) gives the least processing delay compared with the Bartlett, and the rectangular is worst.

R Code Only Below

#41046

require(signal)
N <- 128

r <- rep(1,N)
h <- hamming(N)
b <- bartlett(N)

r <- r/sqrt(sum(r*r))
h <- h/sqrt(sum(h*h))
b <- b/sqrt(sum(b*b))

#r <- r/sum(r)
#h <- h/sum(h)
#b <- b/sum(b)

rd <- rep(0,N)
hd <- rep(0,N)
bd <- rep(0,N)

for (idx in seq(1,N))
{
  if (idx != 1)
  {
    rd[idx] <- rd[idx-1] + r[idx]*r[idx]
    hd[idx] <- hd[idx-1] + h[idx]*h[idx]
    bd[idx] <- bd[idx-1] + b[idx]*b[idx]
  }
  else
  {
    rd[idx] <- r[idx]*r[idx]
    hd[idx] <- h[idx]*h[idx]
    bd[idx] <- b[idx]*b[idx]
  }
}

plot(rd, col="blue", type ="l", ylab="delay", ylim=c(0,max(c(rd,hd,bd))))
lines(hd, col="red")
lines(bd, col="black")

legend(0, 1.0, c("rectangular", "hamming", "bartlett"), col = c("blue", "red", "black"), lty = c(1, 1, 1))


What is the difference between clustering and quantization? The Lloyd-Max quantizer is a scalar quantizer which can be seen as a special case
of a vector quantizer (VQ) designed with the Linde Buzo
Gray (LBG) algorithm.
In k-means clustering, we are given a set of n data points in d-dimensional space  and an integer k and the problem is to determine a set of k points in $R^d$, called centers, so as to minimize the mean squared distance from each data point to its nearest center. A popular heuristic for k-means clustering is Lloyd's (1982) algorithm.
Confusion :

What is the difference between clustering and quantization. I understand the mechanism of clustering technique which is a unsupervised method of grouping data points whereby we map data points into indices of cluster centers which is closest to it. There are different algorithm to perform clustering  - Is clustering another way of doing quantization? Quantization are of 2 kinds -scalar and vector. k-means algorithm is applied for vector quantization. Again k-means algorithm is also applied in vector clustering. So, is vector quantization = vector clustering?
Thank you.

They are essentially the same. The generalized form of the Lloyd's algorithm is LBG which is used for vector quantization (VQ).You may use k-means for the purpose of vector quantization just like the LBG algorithm. In vector quantization, the input space is partitioned into a set of convex areas. They are also referred to as Voronoi regions or cells. If you look at here it says k-means (which is mostly used in computer science) is actually originated from VQ in signal processing as well as pulse code modulation back in fifties.

If you want to implement VQ, you may directly use k-means algorithm. The original LBG algorithm is presented in the following paper if you are interested:

Y. Linde, A. Buzo, and R. M. Gray, ``An Algorithm for Vector Quantizer Design,'' IEEE Transactions on Communications, pp. 702--710, January 1980.


Frequency response equivalent to characterise transient response (to sinusoidal inputs)? I am interested in characterising the time-dependent response of a system to a sinusoidal input.

The steady state response of the system is also sinusoidal, however I wonder if there is a standard way to represent the information in a manner that is more detailed than the Frequency Response Function to capture the transient information too (as a function of frequency input)?

E.g. The system where input: $s(t) = \sin(\omega t), 0 \le t < \infty $ gives output: $y(t) = (1+e^{-t}) \sin(\omega t), 0 \le t < \infty$.

The frequency response would be unity, and the impulse response could not be used to capture the transient behaviour?

As already pointed out in Hilmar's answer, a linear time-invariant (LTI) system is completely described by its impulse response, or - equivalently - by its frequency response, which is the Fourier transform of its impulse response. I think that there is a misunderstanding concerning the Fourier transform, namely that it can only be used to describe steady-state responses to infinitely long sinusoidal input signals, but this is definitely not the case. I think it's illustrative to look at an example:

Let's assume we have a causal real-valued LTI system with frequency response

$$H(\omega)=A(\omega)+jB(\omega)$$

where $A(\omega)$ and $B(\omega)$ are the real and imaginary parts of $H(\omega)$, respectively. Now let's define an input signal to this system:

$$x(t)=u(t)\sin(\omega_0t)$$

Its Fourier transform is

$$X(\omega)=\frac{\pi}{2j}[\delta(\omega-\omega_0)-\delta(\omega+\omega_0)]
+\frac{\omega_0}{\omega_0^2-\omega^2}$$

Using $A(\omega)$ and $B(\omega)$, the imaginary part of the Fourier transform of the output signal $y(t)$ can be written as

$$\Im\{Y(\omega)\}=-\frac{\pi}{2}[A(\omega_0)\delta(\omega-\omega_0)-A(-\omega_0)\delta(\omega+\omega_0)]-\frac{\omega_0B(\omega)}{\omega^2-\omega_0^2}\tag{1}$$

Since $y(t)=0$ for $t<0$ (because this is true for $x(t)$ and because the system is causal), it is fully described by the imaginary part of its Fourier transform:

$$y(t)=-\frac{2}{\pi}\int_{0}^{\infty}\Im\{Y(\omega)\}\sin\omega t\; d\omega,\quad t>0\tag{2}$$

From (1) and (2) we get

$$y(t) = A(\omega_0)\sin\omega_0t+\frac{2\omega_0}{\pi}\int_{0}^{\infty}\frac{B(\omega)}{\omega^2-\omega_0^2}\sin\omega t\;d\omega,\quad t>0\tag{3}$$

Equation (3) fully describes the response of the LTI system to the causal input signal $x(t)=u(t)\sin\omega_0 t$, including transients. This example is supposed to show that the frequency response (or the impulse response) is sufficient for a complete description of the system's response to arbitrary input signals.

EDIT: Derivation of Equation (2):

Let $Y(\omega)=Y_R(\omega)+jY_I(\omega)$ be the Fourier transform of $y(t)$.
We need the following basic properties of the Fourier transform (see here):

$$y(t)\textrm{ real-valued}\Longrightarrow Y(\omega)=Y^*(-\omega)$$
from which follows
$$Y_R(\omega)=Y_R(-\omega)\textrm{ and }Y_I(\omega)=-Y_I(-\omega)$$

If $y_e(t)=0.5[y(t)+y(-t)]$ is the even part of $y(t)$ and $y_o(t)=0.5[y(t)-y(-t)]$ is the odd part of $y(t)$ then

$$y_e(t)\Longleftrightarrow Y_R(\omega)\\
y_o(t)\Longleftrightarrow jY_I(\omega)$$
These relation can be found in the link above.

If $y(t)$ is real-valued the inverse transform can be written as

$$y(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}Y(\omega)e^{j\omega t}d\omega=
\frac{1}{2\pi}\int_{-\infty}^{\infty}\Re\{Y(\omega)e^{j\omega t}\}d\omega=\\
=\frac{1}{2\pi}\int_{-\infty}^{\infty}[Y_R(\omega)\cos(\omega t)-Y_I(\omega)\sin(\omega t)]d\omega=\\
=\frac{1}{\pi}\int_{0}^{\infty}[Y_R(\omega)\cos(\omega t)-Y_I(\omega)\sin(\omega t)]d\omega$$

where the last equality follows from the fact that the integrand is an even function of $\omega$.
Consequently we have

$$y_e(t)=\frac{1}{\pi}\int_{0}^{\infty}Y_R(\omega)\cos(\omega t)d\omega\\
y_o(t)=-\frac{1}{\pi}\int_{0}^{\infty}Y_I(\omega)\sin(\omega t)d\omega$$

With these preliminaries we can finally state the result using the fact that for causal $y(t)$ the following holds for $t>0$ (because $y(-t)=0$ for $t>0$):
$$y(t)=2y_e(t)=2y_o(t),\quad t>0$$

This means

$$y(t)=\frac{2}{\pi}\int_{0}^{\infty}Y_R(\omega)\cos(\omega t)d\omega=\\
=-\frac{2}{\pi}\int_{0}^{\infty}Y_I(\omega)\sin(\omega t)d\omega,\quad t>0$$

The last equality is the same as Equation (2) above.


Tracking an object in a video with I need to track an object in a video with MATLAB. There is an example frame from the video above. I need to track that car that carries the pendulum. But my problem is there is a bar behind the car that has similar color to the car. So I can not just apply img2bw with certain threshold to find the car.

What can I do to isolate it? Since it is my homework I am not asking for a complete code. I want your advice.

I think I can come up with two suggestion to address this problem.

First (simple)

You can simply use HSV color space, and try to find right thresholds. It seems in HSV colour space, the car is pretty distinguishable from the background but the rail is the problem now.

I=imread ( your image );
J=rgb2hsv(I);
imshow(0.1<J(:,:,2) & J(:,:,2)<0.2 & 0.2<J(:,:,1) & J(:,:,1)<0.54,[]);

Second (a bit more complex)

Since your background is presumably static, you can remove background and the remaining would be your car (and of course the rod). I suggest to use optical flow for this purpose.


Interpretation of the eigenvalues of the inverse Hessian in a KLT tracker I'm a masters student, preparing a seminar in computer vision. Among the topics is the Kanade-Lucas-Tomasi (KLT) tracker, as described in

  J. Shi, C. Tomasi, "Good features to track".  Proceedings CVPR '94.

Here's a web resource that I'm using to understand the KLT tracker. I need some help with the math, as I'm a bit rusty in linear algebra and have no prior experience with computer vision.

In this formula for $\Delta p$ (step 5 in the summary), note the inverse Hessian:

$$\Delta p = H^{-1}\Sigma_x\left[\nabla I \frac{\partial W}{\partial p}\right]^\mathsf{T} \left[T(x) − I(W(x; p))\right]$$

In the article, good features to track are defined as ones where the sum of inverse Hessian matrices have large, similar eigenvalues: $\min(\lambda_1,\lambda_2)>threshold$. I was unable to understand how and where this is derived from, mathematically.

The intuition is that this represents a corner; T get that. What does that have to do with eigenvalues? I expect that if the values of the Hessian are low, there's no change, and it's not a corner. If they're high, it's a corner. Does anyone know how the intuition of cornerness comes into play in the eigenvalues of the inverse Hessian in order to determine $\Delta p$ across iterations of the KLT tracker?

I've been able to find resources claiming that the inverse Hessian correlates to the image covariance matrix. Furthermore, the image covariance indicates the intensity change, and then it makes sense... but I've been unable to find what exactly an image covariance matrix is with respect to an image, and not a vector, or a collection of images.

Also, eigenvalues have meaning in principle component analysis, which is why I get the idea for an image covariance matrix, but I'm not sure how to apply this to the Hessian, as it's usually applied to an image. The Hessian, as far as I understand, is a $2\times 2$ matrix defining the 2nd derivatives for $x$, $y$, and $xy$ at a certain location $(x,y)$.

I would really appreciate help with this, as I've been on it for 3+ days, it's just one small formula and time is running out.

Think of them as 2D smoothness terms.
The smoother the patch, the lower the matrix rank and the closer the matrix is to being singular.

On an straight edge (not a corner), just one eigenvalue will be large.
On a corner both will be large.

Using eigenvalues means that the angle of the edge is not a factor, and at any angle, an edge will give just one large e.v.


How to understand multiple peaks in FFT analysis? In FFT analysis, it is often seen that there are several peak frequencies with integer times of frequencies. For example, in the figure below, there are peaks at $0.2, 0.4, 0.6, 0.8$.

The FFT figure below is transformed from the time history of force coefficient $C_{Y2}=F_y/(\rho D U_m^2/2)$ and displacement $Y_2$ on $Y$ direction of a cylinder immersed in a disturbed fluid and attached to a spring.

How to understand this pattern?
If you are familiar with vibration, how to explain these vibration data?

You are looking at harmonics of the fundamental. There's a lot of literature available to read on the topic.


Designing a 2d low pas filter in fourier domain for an image Firstly this is not a homework question. I am working on a project and my professor gave me this task so I can also learn about filters and fourier domain.

I have a picture that is 2800 x 4500 pixels. I need to make a filter that will only give an arbitrary number of "frequencies".

My thought process is this:
1. Transfer the picture in fourier domain using

F = fftshift(fft2(img));

2: Multiply it by a filter that has ones around the centre and zeroes elsewhere.

3: Transfer the image back into x-y domain

The problem is that I don't know how to make the filter since the picture is non-square. I tried making a circle but that didn't work. How do I make an ellipse shape?

In 1-d it is easy because the filter is a vector with ones around the centre.

I tried looking for answer to this but the pictures in question are always square pictures so I cannot apply the same thing. Also nowhere does it say how it works but rather it is just done.

It's been said so many times before that multiplying the Fourier transform coefficients of an image with a mask of ones and zeros in the frequency domain in order to achieve a short-cut to the ideal brickwall filtering in the time-domain is not a recommended method; and almost never preferred due to unavoidable, and unforeseen, errors it would introduce.

That being said, there are algorithms in computer geometry that describe how to approximate a cirle (or ellipse) on a discrete grid like that of a computer monitor.

The following is a circle algorithm that I remember having taken from wikipedia a decade ago. I have no idea on how it works, but just put it here for your sake.

% MIDPOINT CIRCLE ALGORITHM
Nx=256;          % image pixel size
Ny=256;
I=zeros(Ny,Nx);  % image onto which the circle will be drawn

x0 = 128;     % circle centers x-y
y0 = 128;
r = 50;       % circle radius

% ALGORITHM :
% -----------
f = 1-r;
ddfx = 1;
ddfy = -2*r;

x = 0;
y = r;

I(Ny-y0-r+1,x0) = 1;
I(Ny-y0+r+1,x0) = 1;
I(Ny-y0+1,x0+r) = 1;
I(Ny-y0+1,x0-r) = 1;

while x<y

    if f >= 0
       y = y-1;
       ddfy = ddfy + 2;
       f = f + ddfy;
    end%if

    x = x + 1;

    ddfx = ddfx + 2;
    f = f + ddfx;

    I(Ny-y0-y+1,x0 + x) = 1;
    I(Ny-y0-y+1,x0 - x) = 1;
    I(Ny-y0+y+1,x0 + x) = 1;
    I(Ny-y0+y+1,x0 - x) = 1;
    I(Ny-y0-x+1,x0 + y) = 1;
    I(Ny-y0-x+1,x0 - y) = 1;
    I(Ny-y0+x+1,x0 + y) = 1;
    I(Ny-y0+x+1,x0 - y) = 1;

end%for

figure,imshow(I)

The output will be your mask :

I think you can modify it to work for ellipse and some fill in effects...


How to deal with audio with different length in Audio Classification Currently, I am trying to work with the Dataset UrbanSound8K to try some Audio classification. And I got stuck in the preprocessing step already.

Since the audios are of different lengths, like 4 seconds or 0.3 seconds, I found it impossible to directly pass into the whitening algorithms like PCA even after Feature Extraction, using mel-spectrogram/ MFCC.

So my question is what I can do under such circumstance. I was wondering about zero-padding at the end of the shorter sequence. But it seems not working and not going to yield a nice result.

I saw some people using MFCC and summarizing the MFCCs along the time-axis, like mean, variance, kurtosis, skewness.... I think that would work in this case but I just wonder if there are any other ways to do so.

If the signals were of the same duration but obtained using different sampling rates, so that they have different sample lengths, you can resample the shorter sequence with a higher sampling rate to match that of the longer sequence.


Scipy resample, "fourier method" explanation I am seeking explanation of what "Fourier method" is that is described for resampling in documentation of resample() method. I want to downsample an array. I know what decimation is and how it works, I also know what fourier transform is. I cannot however find "scientific backup" so to say to the resample method in scipy. I saw this topic:

Python's $\tt resample$ vs $\tt resample\_poly$ vs $\tt decimate$

I also searched across some books regarding digital signal processing. Am I missing something obvious? I looked in the implementation of resample and I know it performs fourier transform on signal and then takes half of samples needed from the beginning of transformed array and from the end. Then it inverses fourier transform coming back to time-domain. I cannot find scientific description of this method. In the link above, this was mentioned:

https://en.wikipedia.org/wiki/Trigonometric_interpolation#Relation_with_the_discrete_Fourier_transform

I don't understand how interpolation is connected with downsampling. Could anyone explain and/or link some resources? Scipy docs has this in method description:

Resample x to num samples using Fourier method along the given axis.

The resampled signal starts at the same value as `x` but is sampled
with a spacing of ``len(x) / num * (spacing of x)``.  Because a
Fourier method is used, the signal is assumed to be periodic.

Thanks in advance.

Suppose you have initially a real-valued sequence x of length N. The function is basically doing this:

To upsample, it transforms to the frequency domain and adds N/2 zeros at the end. Then it transforms back to the time-domain.
To downsample, it transforms to the frequency domain and deletes the second and third groups of N/4 elements (which correspond to the half with the highest frequency components). Then it transforms back.

You can try it on your own running step by step the source code for the resample() function.

The theory regarding these methods can be found in this answer.


Is time varying channels same as circulary symmetric complex gaussian? From the literature i understood that if the channel is non time varying then we can consider it as circularly symmetric complex gaussian distributed. So can we consider the same if the channel is time varying ?

From the literature i understood that if the channel is non time varying then we can consider it as circularly symmetric complex gaussian distributed.

That's plain wrong. Some non-time-variant channels are $\sim\mathcal{CN}$, some aren't.

  So can we consider the same if the channel is time varying ?

No.


How to accurately count the time period (hardware)? Input : Square wave which sweeps over time domain. (about 1 MHz)
Output: The time period of high/low value.

Like the timer does, I would like to know the actual time (in sec) when the input goes high or low as well as the proper way to implement it.

For example, use FPGA or Microcontroller to read it via analog input.
OR there is the tool that straight to the point! Please let me know.

If it is not a noisy signal you could use a technique called "zero crossing" to gate a counter. You may use an FPGA to create accumulators to add up times and down times.


How to scale measured QPSK points to determine magnitude error? I am trying to determine the magnitude error between demodulated I,Q values and reference values on a constellation plot for a QPSK signal. The reference points are at (1,0),(0,1),(0,-1), and (-1,0). Of course as the amplitude of the input signal increases, the measured points move further away from the references. My (limited) understanding is that measured points would be near the expected to make the resulting constellation useful - where the proximity is affected by the quality of the signal, but not the amplitude.

Is there a standard scaling method used to shift the measured output to the range of the references for this purpose?

Even though have AGC etc in your system, it is not ensured that the received constellation points will lie exactly at the points that you would expect from the TX constellation.

There can be residual offset, gain and phase rotation (at least). In the following , there's a simple method to estimate and remove offset and gain. Essentially, you calculate the mean of the RX constellation, which is your estimated offset. Then, after offset removal, you measure the energy of the constellation and scale the constellation accordingly, such that the constellation has same energy as tx constellation energy:

N = 10000 # number of TX symbols
QAM = np.array([1+1j, 1-1j, -1+1j, -1-1j])
txEnergy = np.var(QAM)

tx = QAM[np.random.randint(4, size=(N,))]

chanGain = 2 # channel gain (combined of AGC, channel and all other processing)
offset = 0.3-0.5j # channel offset
sigma2 = 0.3  # noise variance

noise = np.sqrt(sigma2/2) * (np.random.randn(N) + 1j*np.random.randn(N))

rx = chanGain * tx + offset + noise  # emulate effects of channel + AGC + ...
# RX is the received signal after demodulation

# estimate mean and remove it from the data
rx_mean = np.mean(rx)
rx_mean_correct = rx - rx_mean

# estimate the energy of the constellation
rx_energy = np.var(rx_mean_correct)

print ("Estimated vs real offset: %s <-> %s" % (rx_mean, offset))
print ("Estimated vs real energy: %s <-> %s" % (rx_energy, chanGain**2*txEnergy))

# scale rx constellation to have same energy as tx constellation
corrected = rx_mean_correct / np.sqrt(rx_energy) * np.sqrt(txEnergy)

# calculate constellation MSE
MSE = np.var(rx_mean_correct - tx)
print ("RX MSE: %f" % MSE)

# plot the stuff
plt.plot(rx.real, rx.imag, 'x', label='Received')
plt.plot(corrected.real, corrected.imag, 'x', label='corrected')
plt.plot(tx.real, tx.imag, 'x', label='tx')
plt.legend();

Program output

Estimated vs real offset: (0.300844347234-0.50994604822j) <-> (0.3-0.5j)
Estimated vs real energy: 8.29630318416 <-> 8.0
RX MSE: 2.297161

In order to remove the phase rotation, you need to have some information about the transmitted symbols. I.e. after the received constellation has been processed as above, you need to compare the RX constellation with the known TX points to get the rotation between them. Then, you rotate the RX constellation back according to the estimated phase rotation.


Machine Learning - any suggestions to solve Python rounding errors? I am working with Python to isolate voiced segments from music using the Jamendo corpus for singing detection. Training a model, I break my audio into frames, and have a label (0,1) for each frame.

Unfortunately, due to rounding errors, my labels fall always 1 or 2 frames short.When extracting MFCC features from ~Train Music/03 - sample.ogg, my output is shaped as (13,3709).

This seems reasonable because the length of that file is 3797322 samples and with a hop_length of 1024, the number of frames works out slightly above 3708.

More details about the MFCC feature extraction in the following code snippet:

def MFCC(audio):
    s = []
    for y in audio:
        mfcc = librosa.feature.mfcc(y= y, sr = 16000, n_mfcc=13, n_fft=2048, hop_length = 1024)
        s.append(mfcc)

    return(s)

To now "mark" each frame with the class label (1,0 for singing, no singing respectively), I am converting from seconds to frame numbers but get a maximum value of 3708 (for the longest duration in the file). This is done as per the following function:

def labelToFrames(textCorpus):
    output = []
    for block in textCorpus:
        block_start = int(float(block[0]) * 16000 / 1024)
        block_end = int(float(block[1]) * 16000 / 1024)
        singing = block[2]
        block_range = np.arange(block_start, block_end, 1)  # Step size is 1 (per frame number)
        for x in block_range:
            ms_start = '{0:.3f}'.format(x * 1024 / 16000)           # Converting back to milliseconds
            ms_end = '{0:.3f}'.format((x * 1024 / 16000) + 0.064)   # Converting back to milliseconds
            add_to_output = [ms_start, ms_end, singing]
            output.append(add_to_output)
    return output

I guess librosa's mfcc function is more complex than my 16000/1024 calculation (?). I have tried using Decimal, math.floor and also math.ceil within my block_start and block_stop variables, but I can't seem to match my audio frame length. Are there any other techniques recommended to avoid this error?

EDIT:

To address some discussion points:

y,sr = librosa.load("~Train Music/03 - sample.ogg", 16000, mono = True) followed by y.shape, it yields (3797322,). d = np.abs(librosa.stft(y, n_fft=2048,hop_length = 1024, win_length = 2048))**2 gives a shape of (1025, 3709). Finally performing MFCC on that gives (13, 3709)

Sticking with 44.1khz, the label data is referenced 10 frames short.

I am not sure what exactly is amiss here (referring to the code) but here is a list of things that might help:

Librosa's mfcc returns a matrix that is n_mfcc x ceil(len(y)/hop_length), so no surprises there.
Divisions between integers can return reals in Python 3 but not Python 2. In any case, if you have a division of 16000/1024 try to enter it as 16000/1024.0 to make sure that what is returned is a real, even if it is going to be rounded in any way later on. In general, work with reals all the way up to the point where you need the rounding and then round at that point.

That's not so important around block_start but it would be around ms_start.

To find the start of a frame containing time offset $to$, you need to $Frn_{start}= \lfloor{\frac{to}{\frac{NFFT-HopLength}{Fs}}}\rfloor$. Where $\lfloor{\cdot}\rfloor$ denotes floor.

Which, yes, it is essentially what you are doing here since int will floor.

Assuming that there is no further mismatch between the hardcoded sampling rates and the actual file that is read from the fisk or a mismatch between the audio file and its text file containing the annotations, I cannot see where is the problem with the labeling.

Can you share a specific example where the frame number is off by 1 or 2?

Hope this helps.

EDIT:

Since the annotation file is a simple text file, can you please collect the files that seem to be having issues and confirm that the annotations are indeed "correct"? To do this you can use Audacity. Audacity accepts simple tab delimited text files with annotation information. The Jamendo corpus includes annotations in space-delimited format and files ending in .lab. The conversion is ultra simple using gnumeric, to import the file and export it as tab delimited. Once you do that, open the MP3 in Audacity and import its labels from the newly exported text file.

In the following screenshot, I have imported the emporte-moi.mp3, converted it to mono and downsampled to 16kHz. This did not seem to have affected the annotations, although certain points, I would have annotated differently because it seems to "cut" when the voice still fades out. But that's a detail probably.

Would be useful to see if there is something going on with the annotations of the files that seem to be "out of sync" with the frame count.


Going from digital to analog filters I've been learning about digital signal processing for some time, and I understand the principles in the digital world, as well as the conversion of continuous filters in the s-domain to digital filters in the z-domain.

I'd love to learn more about the process of designing analog filters and converting this into circuits that I can try and test. Since I cannot brag about my knowledge in electronics, with the exception of a little experience dabbling with Arduino projects, I don't see how an equation in the s-domain can be translated to a circuit. From what little I have found, most sources go the other way, taking a circuit and analysing it in the s-domain, and the equations that are given usually aren't explained how they came about. This makes it somewhat hard to understand the relationship between the math and circuits. As far as I understand, the sources generally have equations readily available for simple circuits, and rely mostly on theory of parallel and serial filters to derive larger circuits — is this correct? If so, is this the be all end all method of creating larger circuits, or what is the "correct" approach?

My goal with this is to 1) gain a better understanding of designing analog filters, and 2) take my own filters and make circuits for say guitar pedals and similar applications.

To be honest, I'd love to take any shortcut methods (minimal reading, maximum hands-on and results for motivation), such that I can get some circuits going right away, and then later take a deeper learning approach later on. Similarly, any "catches" or "gotya's" that I need to know about circuits, such that I don't end up blowing myself up, would be nice to know about ahead of time.

What is your recommended learning route to get started as fast as possible? Do you have any recommended material (articles, books, videos, etc.) for a speedy start?

The main trick is knowing how typical lumped circuit element impedances are represented in the $s$-domain. Recall the current-voltage relationship for a capacitor:

$$
i(t) = C \frac{dv(t)}{dt}
$$

Transforming this to the Laplace domain yields:

$$
I(s) = CsV(s)
$$

Note that voltage and current still have a proportional relationship in this domain, just like the time domain. This allows us to write a (frequency-dependent) $s$-domain impedance for the capacitor. This impedance can be used with Ohm's law and techniques like nodal analysis to calculate what a circuit's Laplace-domain output would be for a given Laplace-domain input signal (e.g. an impulse, whose Laplace transform is unity).

For a capacitor, this impedance is:

$$
Z(s) = \frac{V(s)}{I(s)} = \frac{1}{Cs}
$$

This gives capacitors their known highpass characteristic; for $s=j\omega$, as $\omega$ gets larger, the impedance gets smaller. For the limiting case of $\omega=0$, $Z(0) = \infty$; a capacitor passes no DC current.

We can follow the same procedure for an inductor, starting with its current-voltage relationship:

$$
v(t) = L\frac{di(t)}{dt}
$$

Apply the Laplace transform:

$$
V(s) = LsI(s)
$$

and again, we can calculate the inductor's $s$-domain impedance:

$$
Z(s) = \frac{V(s)}{I(s)} = Ls
$$

As the dual to capacitors, inductors have a lowpass characteristic. For small frequencies ($s \to 0$), $Z(s) \to 0$ as well.

Given these impedances, it becomes very straightforward to apply basic circuit analysis techniques to determine the $s$-domain transfer function of a network containing resistors, capacitors, and inductors.


How do I know quantitatively if the correlation of two time series is significant? I computed the correlation coefficient of two time series of daily observations, x and y, but noticed that the more sampling points I used in the calculation, the lower the value of the correlation coefficient, until it reaches saturation at very long time.

It looks that the correlation coefficient is related to the number of points and it's unable to set a fixed threshold to the correlation coefficient. So, how can I know quantitatively if the correlation is significant?

I believe you are looking for the normalized cross-correlation coefficient (source), which in your case would be computed as follows. First, you de-mean your time series $x(t)$ and $y(t)$ such that

$f(t) = x(t) - \bar{x}$, and
$g(t) = y(t) - \bar{y}$

You can treat your time series as vectors $\mathbf{f}$ and $\mathbf{g}$. Assuming the vectors are of equal length, the formula for the normalized correlation coefficient is:

$R_{\text norm}  = \frac{\mathbf{f}^T \mathbf{g}}{\| \mathbf{f} \| \| \mathbf{g} \|}$

This is bound to have a maximum value of $1$ (meaning the time series are very correlated), as you are essentially multiplying together two unit vectors.

However, note that if you have different length input vectors, then you should not scale the output.


Avoiding edge repeats in a FFT convolution I am trying to use FFT as a fast way to convolve a image, however one of the downsides of this is that the "kernel" image will repeat when reaching the edges of the sample.

Is there a way to avoid this to make it more like a classic convolution, that cuts off at the edges?

Zero-pad the edges of the image, then truncate the excess edge from the result


Relationship between frame rate and sampling frequency in speech signal processing? While I was going through the short time processing of speech signals, I got stuck at the relationship between the frame rate and the sampling frequency.

Is there any relationship between them for an audio signal so that the change in sampling frequency will have an effect in the frame rate? If yes, what is mathematical relationship between them for given values of window length (L), window shift (R), Sampling frequency (Fs) and Frame Rate (FR)?

The sample rate divided by the frame rate will give you the window shift in samples.  The ratio of the window (length - shift) to the window length will give you an overlap factor.  The sample rate divided by the window length in samples will give you the FFT result bin spacing.


finite impulse response filter I am not specialist in DSP,I am dealing with geophysical time-series(electromagenteic sources).After converting to frequency domain,my friend is trying to filter with FIR(45-65Hz) to get rid of AC noise(50Hz).The questionis about the harmonics,how will FIR filtering influence harmonics?

The best way to get rid of 50Hz AC noise is with an Adaptive filter (usually the LMS algorithm will work fine), but you need to be able to take a reference signal of only the AC noise source, normally this is usually possible because if your system is being contaminated with AC it means you have some electrical equipment nearby or an outlet where you can take a sample.

If you do not want adaptive filter, then you can try a notch filter at 50Hz. I believe what your friend was trying to do is put a lot of zeros in the Region of 45-65Hz. It depends on the frequency response of the filter but normally the Harmonics will not be attenuated to the same degre, and normally AC is not very clean especially if you have transformers in the middle and such.

Adaptive filter is the way :)


OFDM intra-block multipath interference We know that OFDM eliminates the inter-block intereference. But there are still multipath interference in a single OFDM block(symbol). How does LTE handle this intra-block interference induced by multipath components so as to successfully demod. the signal ?

By using cyclic prefix and its associated demodulation techniques, LTE OFDM (nearly) eliminates the inter-symbol inteference (ISI).
The intra-block inteference, if I understand well your question, refers to intra-carrier inteference (ICI) when the orthogonality between sub-carriers is not maintained due to Doppler, delay spread too large, synchronisation error, ...

Ignoring equipment default that can be mitigated by using sub-carrier pilots, all remaining errors are expected to be fixed by powerful channel coding. Of course the less errors before entering to channel decoder, the better the final result is.

Please note that LTE sub-carrier spacing 15KHz and 7.5Khz are well-chosen so that it is smaller than the coherence bandwidth of typical power delay profiles. Then each sub-carrier is in flat-fading, no ICI if orthogonality is maintained.

Another thing to note that even with the latest modern channel codes, a very long codeword are required to achieve comfortable performance. Inteleaver helps to reduce the length.

So, to be short, after eliminating ISI, remaining errors are corrected mainly by channel code, with helps of pilots, inteleaver and parameter dimensionning.


Cross Domain Equivalent to Nyquist Sampling Theorem? In attempting to answer this question by @Oliver here: What characterizies 'causality' for a finite FFT? I have considered the minimum requirement to avoid time domain aliasing in the Discrete Fourier Transform, or more generally any application where the frequency domain is sampled. Similar to sampling in time at at least twice the highest frequency to represent the spectrum without the effects of aliasing, I suggest using a time duration that is at least twice the response time of the underlying continuous time signal to represent the continuous time domain signal (in the DFT) without the effects of time aliasing. Or when the time domain process is restricted to known causal processes, the time duration is at least as long as the response time.

This is the equivalent to Nyquist's Sampling Theorem in the frequency domain; ultimately "sampling in frequency" such that the duration of the time domain waveform is greater than twice its response time.

I understand that the same theory would apply, but given that Shannon in his paper provides the Nyquist theorem in the time domain specifically has made me curious if this property may go by other formally named theorems in other domains?

To illustrate this graphically consider the drawing by RBJ below except replace the frequency axis with the time axis.

Public Domain, https://commons.wikimedia.org/w/index.php?curid=1065579

I'd say that this is not only "similar to a cross-domain equivalent to Nyquist's Sampling Theorem", but it simply is the sampling theorem. The sampling theorem does not specify the domains of the signals involved; it is rather a mathematical condition that a function of a continuous variable needs to satisfy such that it is perfectly represented by equidistant samples. It is irrelevant if the independent variable of that continuous function is time, frequency, space or anything else.


DC Spur in DAC output I have a design in which a signal from ADC is filtered and then played by DAC. In FPGA the real signal is converted to quadrature using DDS with 70 MHz frequency. After filtering, the baseband signal is sent to DAC. DDS of DAC is also set to 70 MHz.

The problem is that there is always a spur in output of DAC at 70 MHz. Noise level is almost -93 dBm and this spur is -82 dBm while the main signal level is -15 dBm.

Level of spur does not depend on the signal level, as long as signal passes through filter. My first guess was that this is because of signal truncation. There are two steps of truncation in design. I have used convergent rounding instead of truncation but the spur is still there; decreased by almost 6 dB but still there. then I wondered that this might come from ADC, but changing the frequency of input signal does not change the spur; It's always at 70 MHz. Then I taught that it might be because of DAC LO but when the signal is in rejection band of filter, the Spur disappears.

Do you have any idea what makes this spur and how can I remove it?

Any help would be appreciated.

First question is your DAC AC-coupled? In that case, the AC coupling will remove the DC offset caused by your DDS.

Second question, unless your DAC is perfect, your DAC will also have a DC offset. Do you know what's the typical DC offset level of your DAC? Again, if the DC offset caused by your DDS is less than say 10% of the intrinsic DC offset of your DAC, I don't think it's worth it to fix it.

Third question, is the ratio between the sampling frequency and the signal frequency a whole number? If it's an even number, you should not have any DC offset, but if you have an odd number of samples per period, you could get a DC offset, For example, with a ratio of 5 because you might get two round-ups and three round-downs giving your a negative DC offset.

Here's what you could try:

change the signal frequency so that the ratio between fs and f is not a whole number.
Try to add phase dithering, this could remove your DC "spur"


What is diffuse field equalisation? I am using impulse responses from http://sound.media.mit.edu/resources/KEMAR.html

These are available as either raw or diffuse-field equalized HRTFs.

My question is: Which do I want? Raw, or diffuse-field equalized?

I am using these IRs to project single notes (generated synthetically using a ring resonator) to 12 locations around the listeners head.

π

Depends mostly on your playback headset. If it's diffuse field equalized, use the diffuse field set. For most headsets, this will be the best choice.

The only difference is overall equalization, so you can turn one into the other with a simple filter that gets applied to all HRTFs


Do I need to worry about synchronization in a digital signal processor? I understand that a digital signal processor consist of multiple MAD (muliply add units), which execute in parallel to implement a difference equation, like this one:

$$ y[n]=10y[n-1] + 45y[n-2] + 45y[n-4]+ 67x[n]$$

and this one;

$$ y[n]=10x[n-1] + 45x[n-2] + 45x[n-4]+ 67x[n]$$

My question is do I need to be careful about synchronization as we do when programming a multicore CPU/GPU and FPGA?
If not why?

It's hard to come up with a general answer to such a question, because CPU architectures, especially those specialized for signal-processing applications, are all different. In most cases, you shouldn't have to worry about execution-unit synchronization on processors that support parallel execution of some combinations of instructions. In order to achieve such parallelism, there are often some criteria that must be met, such as:

Inter-instruction dependencies (or data hazards), where the output of one instruction is the input to a subsequent one, can cause pipeline stalls if the previous instruction's latency is too large.
Often, you'll need to ensure that your input data is stored in a memory that is fast enough to support parallel execution. As an example, on TI C5000 architectures, there are two flavors of internal RAM: DARAM (double-access RAM) and SARAM (single-access RAM). In order to fully utilize the dual MAC units on the processor, your input data must be in DARAM, or the memory interface won't be fast enough to feed the execution units.

If you meet these processor-specific constraints, then typically, the processor's pipeline logic is smart enough to know this and automatically forward instructions to the appropriate execution units, so that the process is transparent to you. It might be more involved than this on more complicated architectures: I seem to recall that some TI C6000 DSPs, for instance, had a pipeline with less protection, so by ordering instructions incorrectly, you could get undesired operation (whether that was reduced performance or incorrect answers, I don't remember).


Face Recognition: Simplistic Explanation on PCA Eigenface Algorithm I m working on a project that I have to use eigenface but I have some uncertainty and I dont know how to deal with it. There are some tutorials about it on internet but I can't understand what exactly they mean.

This is what i know:

First of all you have to make image matrix to a vector. So you have to attach next rows one after another to the first row.
Then if we have many images of one person, we add up every certain index of all images and divide it count of those images. so
mean[1]= image1[1]+image2[1]+..../images count

Questions:

After that we have make deviation image. First question is we do it for every row image? So:
divImage1[1]=image1[1]-mean[1];
So we have to make covariance matrix of deviation image; so:
cov(image1)=cov( transpose(devImage1)* devImage);
Next question is: it's a big matrix.

How should I deal with it? I read that you can just calculate subImage of it. Or as it's a symmetric matrix, you just have to calculate befor diagonal elements. but im not sure its the main idea
So I have to calculate eigenvector of this matrix. And this is our eigenface. So we have one eigenface for every image of one person?

Next question is: I don't know when they give us a new image how we should work with these eigenfaces to recognize that is some one?

The step-to-step explanation in Eigenface seems quite clear to me.

A covariance matrix is like an high-dimensional extension of the variance, which is computed by removing the average from your only sample.

Yes, you remove the average face ($ \operatorname{AF}$ from all images, but keep it preciously. Your cov(image1) definition seems weird to me, but suppose you have it.
The eigenvectors of the big matrix are eigenfaces ($ \operatorname{EF}$), put in vector form, you can reshape them to image form. Imagine that as a common trait shared by several persons in your image base.

At this stage, each image $I_k$ can be approximated by a linear combination of $ \operatorname{EF}$s:
$$I_k =  \operatorname{AF}+a_k \operatorname{EF}_1+b_k \operatorname{EF}_2+c_k \operatorname{EF}_3\ldots\,.$$
The idea is that the different pictures of the same person will have almost the same $a_k,b_k,c_k$, because their share the same traits. But a different person would have quite different coefficients. So comparing the vectors of coefficients may help distinguish different persons. In your case:

You have computed eigenfaces. A new picture arrives. If you compute its coefficients on the $\operatorname{EF}$s (called a projection), if they are close to the coefficients of another person, in fact the closest, then, maybe, the new picture is a picture of the guy.

Additional links:

Eigenfaces, for Facial Recognition
Eigenfaces
Face recognition 101: Eigenfaces


Is it possible to discretely sample the function I had a few questions on sampling(I'm quite new witht his), I tried to answer them, I think that I did the first one correct , but not sure about the 2 other:
.
given the next functions,Is it possible to discretely sample the function? if so what is the maximal allowed distance between the samples? otherwise explain why?

$$f(x)=\sin(\alpha x)$$

For this one I said it is possible using Nyquist theorem, assumin $\alpha =2\pi f$ and $T=\frac{1}{f}=\frac{2\pi}{\omega} $ then the allowed distance is $\frac{T}{2}=\frac{\pi}{\omega}$

(hopefully i got this allright).

now from the second question I'm not sure.

$$f(x) =
\left \{
\begin{array}{cl}
1,  & \text{$-1\leq x \leq 1$} \\
0, & \text{else}
\end{array}
\right .
$$

for this function I'm pretty sure the answer is it is not possible, since it is a straight line, so I'm not sure.. (please help me with this one).

$$f(x)= \mbox{convolution between function of question 1 , function of question 2.}$$

In this case I'm guessing it is possible since its the same as sampling first question function (only that it is sliced).

yes, with spacing between samples closer than $\frac{\pi}{\alpha}$.
no.  it's not sufficiently "bandlimited"
yes, the convolution of a sinusoid with anything is a sinusoid of the same frequency.  so it's "yes" for the same reason as 1.


Difficulty with a Fourier Transform What would be the best way to take the Fourier transform of
$$
f(t)\cdot \cos\big(\pi(t-1)\big)
$$
I'm aware that when you take the Fourier Transform of $\cos(kt)$ you get two impulse at the location of $k$. And that the $(t-1)$ would create a complex exponential. But you can only apply the time-shift when both functions are time shifted (at least to my knowledge at least). Really the $\pi$ and the non-time shifted $f(t)$ are really throwing me off. Thanks!

It is easy to see that using the trigonometric identity in Equation $(1)$ below
$$
\cos(\theta_1 \pm \theta_2) = \cos(\theta_1)\cos(\theta_2) \mp \sin(\theta_1)\sin(\theta_2) \tag{1}
$$
We have
\begin{align}
f(t)\cos\big(\pi(t - 1)\big) &= f(t)\overbrace{\cos(\pi t - \pi)}^{\text{use Equation}\ (1)}\\
&\equiv -f(t)\cos(\pi t)\tag{2}
\end{align}
Which can be rewritten as
$$
-f(t)\cos(\pi t) = -f(t)\cos\left(2\pi \frac t2\right) = -f(t)\cos(2\pi f_c t)\quad\text{where}\quad f_c = \frac 12
$$
And we know the Fourier cosine modulation frequency-shift property in equation $(3)$:
$$
\mathcal F\big\{x(t)\cos(2\pi f_c t)\big\} = \frac 12\big[X\left(f - f_c\right) + X\left(f + f_c\right)\big]\tag{3}
$$
You can now easily use Equation $(2)$ and $(3)$ to find the solution.


Dummy DSP student request for IDTFT Calculation and Integration Process We have one DSP course in biomedical engineering, so i have been for 10 years away from mathematical calculation and working in this years by programing by c and ..., so one of our exercise is to calculating the IDTFT of some signals like this from this eBook "Digital Signal Processing Using MATLAB 3rd Edition":

  Using the definition of the inverse DTFT in (3.2), determine the
  sequences corresponding to the following DTFTs:

so the manual book solved this question by summery like this:

So, i don't know how it reached this part:

could you give me some ...?
It is possible to do this calculation by some software and seeing the formulation changing mode from base question to the answer like those shown above?

Thanks for your attention.

This is using Euler's Identity for a sine which should get you past the point where you are stuck:

$$ \sin(\phi) = \frac{e^{j\phi}-e^{-j\phi}}{2j}$$


Analytical Calculation of the Impulse Response of an Arbitrary Filter This is a rather general question.  Suppose I have an arbitrary filter, $F$, defined in the Laplace domain;
$$
F=k\frac{\prod_i (s+z_i)}{\prod_j(s+p_j)}.
$$
How do I calculate the impulse response of this filter analytically?  As a specific simple example, consider the first order low pass filter
$$
F_0=\frac{1}{s+p}.
$$
To be clear, what I want to get from this transformation is a function which describes the delta function response from $t=0$ to $t\rightarrow\infty$ where the delta function is applied at $t=0$.

Your expression is the Laplace transform of the exponential function $e^{-pt}$. The Laplace transform of the dirac unit is $1$. So the inverse Laplace transform is exactly $e^{-pt}$. Maybe the fact that this function is noncausal adds a little complexity here.


Nyquist frequency and DC Im studying the DSP book by Steven W. Smith. On page 41 he covers the Nyquist frequency. He makes up an example by writing the following.

  ... consider an analog signal composed of frequencies between DC and 3kHz.

What does that mean? My guess is that in this context DC is equivalent to zero. Why does he use DC here?

Sorry if this is totally obvious.

You are right, DC means zero frequency. The term probably originates from electrical engineers' use of the terms DC and AC.


basic bandpass filter I am attempting to implement a basic bandpass filter but the center frequency of my filter seems to be irrelevant as I can change it all I want but it has no effect on the filter.
where am I going wrong?

I have given my calculations here in the hope that someone will be able to point out where I'm making a mistake.

Filter specs:
supression of DC and Fs/2 (zeroes at +1 and -1)
center frequency at pi/4
bandwith of pi/16

2(1-R)=pi/2
1-R=pi/32
R=(pi/32)-1
R=-0.9018252296

 K*(z-1)*(z+1)
---------------
z-R*cos(pi/4)+R

     K*(z^2 -1)
-------------------
z^2-R^2*cos(pi/4)+R

all leading to a difference equation of

Y[n]=R^2*cos(pi/4)*Y[n-1]-R*Y[n-2]+K*(X[n]-X[n-2])

Your denominator is wrong. It should be on this form

$$
(1-z_0z^{-1})(1-z_0^*z^{-1}) = 1-2R\cos(\phi)z^{-1}+R^2z^{-2}
$$
where $z_0 = Re^{j\phi}$.

If your goal is to block DC and create an anti-aliasing filter in a simple way, I would recommend not to use complex poles. Instead use the procedure described in detail here: http://www.dsprelated.com/showmessage/172787/1.php.

It even provides a C implementation.


Calculating RMS/Crest Factor for a Stereo Signal I'm doing some MIR related work with stereo audio and am a little unsure on how to proceed with some basic tasks.  All my academic work prior to this has involved forcing audio to "mono" to perform analysis, and I'm trying to avoid that as much as possible here.

So, for example, I have a crest factor function in python that looks like this for an array of audio data (either mono, left, or right channel):

def get_cf(data, win_size):
    """
    data: audio array in mono, left, or right channel only
    win_size: size in samples for the block analysis (created in calc_crest_factor)

    calc_crest_factor passes mono style data to this function to get the crest factor.

    returns: the crest factor for each window"""
    # Buffer the signal matrix-style (input, block-size, hop-size)
    data_matrix = librosa.util.frame(data, win_size, win_size)

    peaks = np.amax(np.absolute(data_matrix), axis=0)

    # Get the mean-square over each window
    RMS = np.sqrt(np.mean(np.square(data_matrix), axis=0))

    # Get crest factor for each window
    return np.divide(peaks, RMS)

In another function, I call the function above like so for stereo audio:

 if len(data) == 2:
    crest_factor_l = get_cf(data[0,:], win_size)
    crest_factor_r = get_cf(data[1,:], win_size)

And here I now have a crest factor for the right channel and one for the left channel, but am unsure on how to proceed on getting a final "crest factor" value for the audio.  Do I just choose the greater of the two values for each window?  Or is it more correct to keep both values?

This same question goes for RMS too.  I know in my pro audio applications, stereo audio gets passed to some RMS calculator and a single value gets pumped out.  How do the two separate channels with two separate values for either Crest Factor or RMS become 1 value while honoring the "stereoness" of the data?

Conceptually you can think of it as creating one mono signal by concatenating the left and the right channel signals, and then computing the RMS and the Crest factor of that mono signal. With $N$ the number of samples in the current frame, and $x_l[n]$ and $x_r[n]$ the left and right channel signals, respectively, you get the following formulas:

$$\begin{align}x_{\text{rms}}&=\sqrt{\frac{1}{2N}\left[\sum_{n=0}^{N-1}x^2_l[n]+\sum_{n=0}^{N-1}x^2_r[n]\right]}=\sqrt{\frac12\left(x^2_{l,\text{rms}}+x^2_{r,\text{rms}}\right)}\\
C&=\frac{\max\{|x_l[n]|,|x_r[n]|\}}{x_{\text{rms}}}\end{align}$$

where $x_{l,\text{rms}}$ and $x_{r,\text{rms}}$ are the RMS values of the left and right channel, respectively.


MRI image reconstruction - reduced scan pecentage I have an MRI kspace from a siemens 1.5 T scanner with reduced scan percentage ( sampled truncation by not acquiring the most peripheral lines in kspace). The kspace plot in matlab look like below.

I can do zero filling in this kspace to the required resolution and reconstruct images.

Other than zero filling, is there another method which can estimate missing lines and reconstruction images.

I would say no, there isn't. This is unmeasured information and you cannot compute it without additional information (ie some additionally sampled data, and maybe some compressed sensing algorithm). However, truncating an image like that in order to save measurement time is quite unusual. You would be better off to sample one side full and skip twice the number of lines on one side of $k$-space. You could then use additional information (in this case: the object is real and should be shimmed well, ie you know that the image data should be real). The $k$-space of a real image is Hermitian, and you could use the sampled information of the one side outer part to use the Hermitian symmetry to fill the missing lines. Different algorithms for this exist.


What does photon noise mean in digital image processing or computer vision? From this link:

  Photon noise, also known as Poisson noise, is a basic form of uncertainty associated with the measurement of light, inherent to the quantized nature of light and the independence of photon detections.

I couldn’t understand,

  … and the independence of photon detections.

What did the author actually mean by that?

Ultimately at a fundamental level, an image signal as captured by a digital sensor is obtained by integrating the energy transmitted by the photons of the light wave.

If you are imaging a bright scene, there will be a large enough number of photons hitting the detector that this effect will not be detectable, but if the sensor is operating in a dark environment, individual photons may be detected and counted, and the number of photons that have hit a cell of the detector as a function of time is modeled as a Poisson random process.


ASK demodulation using DSP I need to implement an ASK (ampliltude shift keying) demodulation using a DSP ( a dsPIC ) and programming in C. Especifically talking, I need to demodulate a signal from a car's key transponder. I already have a circuit of antenna that receive the signal and I already can sample it, but I still have a doubt of how could be the best way, faster and reliable, to filter the signal.
I was using a FIR filter using Kaiser windowing and 13 taps, but it demanded a lot of processing time and I still have some errors. Is there a better way to do it ?

With my best regards,
Daniel.

In general, a receiver would consist of these blocks:

An analog band-pass filter centered around the carrier frequency.
A low-noise amplifier to bring the signal above the receiver's sensitivity.
Optionally, an analog intermediate frequency stage.
A sampler. You can do bandpass sampling if the carrier frequency is too large for the ADC.
In a DSP, implement:

Final downconversion if needed.
Downsampling to reduce the sampling rate
Carrier, frame and symbol synchronization
Matched filtering to recover the bits with optimum BER (in the case of ASK with rectangular pulse shaping, the matched filter is just an integrate & dump filter)

I highly recommend that you get and read Telecommunications Breakdown, by Johnson and Sethares, which is available for free (as a draft) from the authors' website here.


Can I use FFT to interpret accelerometer gestures? I have the need to detect two different gestures that occur when looking at accelerometer data. Here's a run down (as brief as I can make it):

Lets say an iPhone is being oscillated back and forth while remaining face-up. The  user can either do one oscillation (once forward or backward, gesture one) or a continuous oscillation for any amount of time (gesture two).

Currently my code maintains a list of the device's accelerometer data (y-axis) from the last 50 frames . Each frame this data is put through an FFT algorithm (this one http://goo.gl/yi3mn), and then I try to interpret the frequency domain given. I have noticed a strong relationship between the speed of the oscillations, and the power of the mid to low range of the frequency space.

The problem is that I need to detect (in real-time as the movement animates something on screen) whether the just finished oscillation is a singular one, or one that continues into another oscillation in the opposite direction. Here I should point out that accelerometers work with pressure plates. When an oscillation stops the input data will show the opposite value from when the oscillation was being performed. This makes it difficult to see any obvious difference between a single forwards movement and a double forwards then backwards movement (at the end of the each first oscillation).

Can anyone suggest how I might make use of an FFT to decipher an accurate magnitude value (or something else useful) that I can use to tell the difference between the gestures described above?

I have recorded the data I am using to text files. It is the raw data (g-force values) along the y-axis of the accelerometer input.

http://pastebin.ca/2108123 shows data for 2 singular oscillations (I hold the device still at the start, end, and between the two oscillations).

NOTE: new raw data with 20 oscillations has been uploaded, but not yet plotted.
http://pastebin.ca/2108387 shows data for 20 continuous oscillations (I hold the device still at the start and end).

Looking at the data I can't see why you'd want to use FFT for this purpose (I may be wrong though). Lets say you want to detect single or multiple movements based on this data set (which you've attached) the methods I'd be looking at would be

Something on the lines of rate of change of absolute values. First derivative would do nicely I guess.
Wavelet transform maybe? I can't elaborate more on its utility until I know whether you want to know the the position of the peaks or not.

If the only thing you want to detect is single vs multiple oscillation then FFT is, in my opinion, not a very good choice for analysis.


FIR Impulse response & transfer function I am working though a dsp past paper and I have come across the questions "find the impulse response in the time domain" & "find the transfer function in the time domain" I know its really simple but the "time domain" part of the question catches me out, I know how to find the transfer function in the Z domain, its not for homework, its a past paper I am looking through.

How would I find the impulse response in the time domain?

AND

How would I find the transfer function in the time domain?

I have already tried YouTube and google and a few books but haven't found anything that is specific to those questions. The question is attached in the photo.

I agree with Matt's comment, however i suspect that whoever gave you the assignment wants you to write down an equation describing the transfer behaviour in time domain for an arbitrary input.

When $x[n] = \delta[n]$, the impulse response of the given system would be

$$h[n] = -\delta[n] + 2\delta[n-1] - 2\delta[n-2] \ .$$

You can find it in the output $y[n]$ by setting the input $x[n]$ to

$$ x[n] = \delta[n] \triangleq \begin{cases}
1 \qquad \text{if } n = 0 \\
0 \qquad \text{if } n \ne 0 \\
\end{cases}$$

Then the output at the end of the diagram is $h[0]=-1$. The $z^{-1}$ elements act as delays of one sample, so the next outputs are $h[1]=2$ and $h[2]=-2$, respectively.

Now to find out what the equation does to general inputs in time domain, write the same thing down in terms of $x[n]$ instead of $x[n]=\delta[n]$.

For any given $x[n]$, the resulting $y[n]$ is produced by the current input times $-1$ and the prior two inputs multiplied with their factors $2$ and $-2$. You can then write down the general solution as

$$y[n] = -x[n] + 2x[n-1] -2x[n-2] \ .$$

Keep in mind however that the term "time domain transfer function" could be misunderstood by some, which is probably the reason your search did not yield any results.


Counting shape changing pattern in signal I need to count the number of 'patterns per second' in a signal. The problem is that the signal changes 'shape' a lot and unpredictably. It is still easy for a person to see the the pattern, but I don't know how to do this with code.

EDIT For extra clarity: say the first sample is 1 second. There are 17 copies of the pattern in this sample. So I want to compute "17 (patterns per second)".
Also, it might be worth noting the samples below are from a continuous measurement, in this order.

Any ideas how to tackle this?

(It needs to count the pattern in realtime.)

Here are some samples of the signal. (Zoom level in both directions is identical for all images.)

You can trying using a musical pitch detection/estimation algorithm, which may be more robust in some cases than just a simple autocorrelation (especially in one of your examples that includes a slightly inharmonic overtone).  From the frequency or period of the estimated "pitch" you can compute the number of pattern repetitions per second.

A bunch of pitch estimation methods are mentioned in my answer to this question: https://stackoverflow.com/questions/7181630/fft-on-iphone-to-ignore-background-noise-and-find-lower-pitches


Why the power of periodic signal is square? I have DSP in my academics and while going through the video lectures i am stuck regarding power of periodic signal.

Its mentioned as

$$P_{\text{avg}}= \frac{1}{N}\sum_{n=0}^{N-1}|x(n)|^2$$

my doubt is why we have to mention $|x(n)|$ as square ?

Note: I am having format issues please edit my Question and put it in format. thank you.

Power is proportional to the squared amplitude of a signal. This "power" concept comes from thinking of the signal $x[n]$ as representative of a voltage or current waveform. In that model, the power dissipated by applying $x[n]$ to a one-ohm resistor is the value of $|x[n]|^2$ averaged over one period.


How can I generate a sine wave with time varying frequency that is continuous, in C? How can I generate a sine wave with time varying frequency that is continuous? How can I resolve the following problem? I want a continuous graph.

I'm generating this simply like this:

for(int i = 0; i < pcm_buffer_size - 1; i += 2) {
    float sample = gain * sin(((float) t * M_PI * 2 * 500) / (float) (sample_rate));

    printf("%f\n", sample);

    t++;
}

for(int i = 0; i < pcm_buffer_size - 1; i += 2) {
    float sample = gain * sin(((float) t * M_PI * 2 * 1000) / (float) (sample_rate));

    printf("%f\n", sample);

    t++;
}

Where 500 and 1000 is the frequency.

Well lets go lol

@AnonSubmitter85 give to you a nice answer, but let me show my  way to do it in matlab, and this maybe can be very easy to port to C:

First I'm creating 256 samples in 500hz sampled at 44100hz take a look how I accumulate the phase and in the end of the first loop I put the phase between the interval 0 and 2pi...

Nice now lets go to the second loop to create more 256 samples of 1000Hz and to do it continuous I'm using the last phase :-), here is the code:

m_phase=0;

signal=[];

f=500;
fs=44100

phaseInc = 2*pi*f/fs;

for i=1:256
   signal(i) = sin(m_phase);
   m_phase = m_phase + phaseInc;
 end

%place the phaser between the 0 and 2pi range
m_phase = mod(m_phase, 2*pi);

f=1000
phaseInc = 2*pi*f/fs;

for i=257:256*2
   signal(i) = sin(m_phase);
   m_phase = m_phase + phaseInc;
end

%place the phaser between the 0 and 2pi range
m_phase = mod(m_phase, 2*pi);

here is the plot of the above code::


What should be the cut off frequency of a Filter operating on an OFDM signal? I have a time domain OFDM signal and to this I need to apply different linear filter however I am not sure on the cut off frequency. I am simulating for multiple symbols where each symbol has 1200 central subcarriers in a 2048 FFT/IFFT bin. The samp freq is 30.72MHz. While designing filters in MATLAB using commands like "fir1", I am not sure how to set cut off freq.

You don't need that filter. Your previous questions still indicate you shouldn't be using a filter. OFDM is a filterbank in itself, and applying a filter beforehand will increase delay spread, and thus you get ISI¹.

Sure, with a filter, you'll save on the noise power you'd get in carrier sidelobes, but that should typically be lower than what you lose in ISI.

So, no, you don't apply a different linear filter. You simply ignore the FFT bins that don't contain signal.

¹ you'll get ISI under the assumption that your OFDM system is designed with the minimum amount of subcarriers needed to compensate the maximum expectable delay spread on the channel. That makes a lot of sense because using significantly more subcarriers has no advantage, but makes your system depend on a longer coherency time. If you convolve the channel output with your receive filters, you get a "convolved channel-filter", with a length that is thus longer than your channel, and thus worse than your OFDM system can deal with.


filter coefficients to know high pass and low pass filter Looked on few questions already there in dsp.stackexchange and couldn't get the direct answer. Basically I have a kernel of a gaussian separable filter.

float w[N*2+1] = {-0.00467444956f, -0.0531099625f, 0.0152470656f, 0.300147355f,  0.484779954f, 0.300147355f,  0.0152470656f, -0.0531099625f, -0.00467444956f};

float w2[N*2+1] = {-0.00899997167f, 0.0197822638f, -0.0300691985f, 0.0374629796f, 0.963647842f, 0.0374629796f, -0.0300691985f, 0.0197822638f, -0.00899997167f};

The coefficients in first first increasing and then starts decreasing.The middle point is strongest and it is close to 0.5.

In 2nd the coefficients are of alternate sign with middle being storngest close to 1.

What properties these coefficients signify?
How can I relate them to high pass or low pass filter?

I am using the above kernels in the smoothening the image using Gaussian separable filter.

I think the filters may be sum of a lowpass and a hightpass filter. They may use to extract some high frequency and low frequency components of a signal.
For example to extract the edges of an image and smooth then simultaneously.