Nyquist Frquency isn't what it is supposed to be I am new here and also am not very knowledgeable about DSP so this might be dumb and easy. I am aware of the fact that in order to reconstruct a signal, I need to sample it at a frequency that is more than twice its highest frquency.

In the example below, I use N=50 samples and expect the reconstruction to fail if I use freq = 25 or so (perhaps due to numerical round-off). But I get that the reconstruction breaks around freq = 12.5. I feel I am missing something fundamental.

Here are plots and the code to generate them:

from numpy.fft import rfft as rfft
from numpy.fft import irfft as irfft
import numpy as np
import matplotlib.pyplot as plt
import math

def interpolant( f_hat, pts ):
    '''
    This is the trigonometric polynomial interpolating
    the signal that is encoded in the values of f_hat.
    it is evaluated at the points pts

    The commented part is the true calculation to be done,
    the uncommented just does it more accurately.
    '''
    N = len( f_hat )
    f = f_hat[0]
    for k in range( 1, N/2 ):
        f = f + 2 * f_hat[k].real * np.cos( k * pts )
        f = f - 2 * f_hat[k].imag * np.sin( k * pts )
        #f = f + f_hat[k]             * np.exp( 1j * k * pts )
        #f = f + f_hat[k].conjugate() * np.exp(-1j * k * pts )
    f = f + f_hat[N/2] * np.cos( N/2 * pts )
    return f / N / 2

def f( pts, freq ):
    '''
    this is "the conitinuous" signal, evaluated (sampled)
    at pts
    '''
    return  np.sin( freq * pts )

# Number of sampled points
N = 50

# Where we sample the signal
pts = np.linspace(0, 2 * np.pi, num = N, endpoint = False )

# A much finer grid, used solely to display results
oversampled = np.linspace( 0, 2 * np.pi, 20 * N, endpoint = False )

# The grid we use for the plotting
grid = oversampled

# The frequncy of the sine wave below
freq = 7

# Do the FFT on the sampled signal
f_hat = rfft( f( pts, freq ) )

plt.plot( grid, interpolant( f_hat, grid )  , color = "g" )
plt.plot( grid, f( grid, freq )      , color = "r" )
title1 = str(N) + " samples. Signal frequency is " + str(freq)+" \n"
title2 = "Red is true, green is interpolant. Reconstruction succeeds"
plt.title( title1 + title2 )
plt.show()

# The frequncy of the sine wave below
freq = 12.9

# Do the FFT on the sampled signal
f_hat = rfft( f( pts, freq ) )

plt.plot( grid, interpolant( f_hat, grid )  , color = "g" )
plt.plot( grid, f( grid, freq )      , color = "r" )
title1 = str(N) + " samples. Signal frequency is " + str(freq)+" \n"
title2 = "Red is true, green interpolant. Reconstruction fails"
plt.title( title1 + title2 )
plt.show()

The closer the sampling frequency gets to half the (highest spectral) frequency of the signal, the longer one has to sample to get a sufficient bounded fit (for some bound) of a reconstruction using those finite number of samples.  This required sampling time extends to infinity as one approaches Fs/2.

So, you simply didn't sample for long enough in your experiment.  Try 5 thousand and 5 million samples for sampling a 12.9 Hz sinusoid at 26 sps, and you will see a closer fitting looking plot over most of an X axis of 5 thousand or better million reconstructed sample points.  Then try a sinusoid even closer to Fs/2 and you will find you need to sample even longer to get the same visual quality of fit.


difference between GMSK (Gaussian Minimum shift keying) and GFSK (Gaussain frequency shift keying) I am unable to differentiate GMSK and GFSK.
How do the modulation techniques differ in Tx?

How do their receiving techniques differ..?

With GMSK, the frequency modulation index is 0.5. With GFSK, the frequency modulation index is larger than 0.5 (the signaling tones are farther apart).
The receiving techniques used for GMSK as an FSK work for GFSK.


How to compute the channel gain from the path loss model of a wireless channel I did my research but I did not get a satisfying answer.
We have a transmitter and a receiver over a wireless communication channel. The channel's path loss is modeled as
$ 128 + 37 \log_{10} (d)$, with a shadowing standard deviation of $10$ dBs.

I am assuming that this path loss equation is expressed with dBs (nothing explicitly says that).

I have read that the combined path loss and the shadowing model is just the addition of the path loss model and a normally distributed random variable with mean $0$ and std. $10$ dBs.

Is this addition of path loss and shadowing is the same as the channel gain ?
If not, how to calculate the channel gain $\vert h \vert$ from the path loss model, combined with the shadowing ?

To answer the literal questions you asked:

yes
see above

To Frame Challenge your question and give you what you were probably looking for:

The path loss equation you cited:
$$
128+37\log_{10}(d) ,
$$
is equivalent to:
$$
\frac{P_r}{P_t} = \frac{10^{-12.8}}{d^{3.7}},
$$
where $P_r$ is the received power in Watts, and $P_t$ is the transmit power in Watts.  From this equation it can be seen that the path loss exponent is 3.7.

Under log-normal shadowing, the random shadowing component is modeled as a zero-mean Gaussian random variable (in dB) added to the path loss.


How to remove spike from the data series? I am getting load cell reading from my adc ADS1231 every 100 ms or say at 10Hz i.e. 10 samples per seconds. I have obtained the readings as shown in the image.
As you can see i am getting spike on sample number 8 which I would like to eliminate using digital filtering method. Any advice would be of great help.

I am continuously getting the adc values and the older data will be the first one to go out  the accumulator. In below image i have shown accumulator values from 1 to 10. Latest value will be at position 10 and oldest will be as position 1.

I would put this in a comment but I don't have enough rep points, as just started on this.

I would recommend reading this information before designing with sigma-delta ADCs
http://www.numerix-dsp.com/appsnotes/APR8-sigma-delta.pdf

It may be that you need a basic analog filter before the ADC input to prevent interfrence noise aliasing back in to your signal (assuming the load on your cell is constant?).


Finding symmetric regions/patterns in image I have a set of images representing the mean curvature of a human back surface.

What I want to do is to "scan" the image for points which have similar, reflected "counterparts" in some other part of the image (most likely symmetrical to the midline, but not necessarily since there might be deformities). Some image-stitching techniques use this to "auto-detect" similar points between images, but I want to detect them for both sides of the same image.

The ultimate goal is to find a continuous, most probably curved, longitudinal line that adaptively divides the back in symmetric "halves".

A sample image is placed below. Notice that not all the regions are symmetrical (specifically, just above the center of the image, the red vertical "strip" deviates to the right). That region should receive a bad score, or whatever, but then the local symmetry would be defined from symmetrical points placed farther away. In any case, I'll have to adapt any algorythm to my application domain, but what I'm after is som correlation/convolution/pattern-matching strategy, there must already be something around, I think.

(EDIT: there are more images below, and some more explanation)

EDIT: as requested, I'll include more typical images, either well-behaved and problematic. But instead of colormapped images, they are grayscale ones, so that color relates directly to data magnitude, which didn't happened with the colored image (provided just for communication). Although the gray images seem to lack contrast compared to the colored one, the data gradients are there and can be brought up with some adaptive contrast if desired.

1) Image of a very symmetrical subject:

2) Image of the same subject at a different moment. Although there are more "features" (more gradients), it doesn't "feel" so symmetric as before:

3) A thin young subject, with convexities (bony protrusions, denoted by lighter regions) at the midline instead of the more common concave midline:

4) A young person with spinal deviation confirmed by X-Ray (note the asymmetries):

5) The typical "tilted" subject (although mostly symmetrical around the curved midline, and as such not properly "deformed"):

Any help is much welcome!

As I've said in the comments, medical image registration is a topic with lots of research available, and I'm not an expert. From what I've read, the basic idea commonly used is to define a mapping between two images (in your case an image and its mirror image), then define energy terms for smoothness and for image similarity if the mapping is applied, and finally optimize this mapping using standard (or sometimes application-specific) optimization techniques.

I've hacked together a quick algorithm in Mathematica to demonstrate this. This is not an algorithm you should use in a medical application, only a demonstration of the basic ideas.

First, I load your image, mirror it and split these images into little blocks:

src = ColorConvert[Import["http://i.stack.imgur.com/jf709.jpg"],
   "Grayscale"];
mirror = ImageReflect[src, Left -> Right];
blockSize = 30;
partsS = ImagePartition[src, {blockSize, blockSize}];
partsM = ImagePartition[mirror, {blockSize, blockSize}];
GraphicsGrid[partsS]

Normally, we would do an approximate rigid registration (using e.g. keypoints or image moments), but your image is almost centered, so I'll skip this.

If we look at one block and it's mirror-image counterpart:

{partsS[[6, 10]], partsM[[6, 10]]}

We can see that they're similar, but shifted. The amount and direction of shift is what we're trying to find out.

To quantify the match similarity, I can use the squared euclidean distance:

ListPlot3D[
  ImageData[
   ImageCorrelate[partsM[[6, 10]], partsS[[6, 10]],
    SquaredEuclideanDistance]]]

sadly, using this data is the optimization directly was harder than I thought, so I used a 2nd order approximation instead:

fitTerms = {1, x, x^2, y, y^2, x*y};

fit = Fit[
   Flatten[MapIndexed[{#2[[1]] - blockSize/2, #2[[2]] -
        blockSize/2, #1} &,
     ImageData[
      ImageCorrelate[partsM[[6, 10]], partsS[[6, 10]],
       SquaredEuclideanDistance]], {2}], 1], fitTerms, {x, y}];

Plot3D[fit, {x, -25, 25}, {y, -25, 25}]

The function is not the same as the actual correlation function, but it's close enough for a first step. Let's calculate this for every pair of blocks:

distancesFit = MapThread[
   Function[{part, template},
    Fit[Flatten[
      MapIndexed[{#2[[2]] - blockSize/2, #2[[1]] - blockSize/2, #1} &,
        ImageData[
        ImageCorrelate[part, template,
         SquaredEuclideanDistance]], {2}], 1],
     fitTerms, {x, y}]], {partsM, partsS}, 2];

This gives us our first energy term for the optimization:

variablesX = Array[dx, Dimensions[partsS]];
variablesY = Array[dy, Dimensions[partsS]];

matchEnergyFit =
  Total[MapThread[#1 /. {x -> #2, y -> #3} &, {distancesFit,
     variablesX, variablesY}, 2], 3];

variablesX/Y contains the offsets for each block, and matchEnergyFit approximates the squared euclidean difference between the original image and the mirrored image with the offsets applied.

Optimizing this energy alone would give poor results (if it converged at all). We also want the offsets to be smooth, where the block similarity tells as nothing about the offset (e.g. along a straight line or in the white background).

So we set up a second energy term for smoothness:

smoothnessEnergy = Total[Flatten[
    {
     Table[
      variablesX[[i, j - 1]] - 2 variablesX[[i, j]] +
       variablesX[[i, j + 1]], {i, 1, Length[partsS]}, {j, 2,
       Length[partsS[[1]]] - 1}],
     Table[
      variablesX[[i - 1, j]] - 2 variablesX[[i, j]] +
       variablesX[[i + 1, j]], {i, 2, Length[partsS] - 1}, {j, 1,
       Length[partsS[[1]]]}],
     Table[
      variablesY[[i, j - 1]] - 2 variablesY[[i, j]] +
       variablesY[[i, j + 1]], {i, 1, Length[partsS]}, {j, 2,
       Length[partsS[[1]]] - 1}],
     Table[
      variablesY[[i - 1, j]] - 2 variablesY[[i, j]] +
       variablesY[[i + 1, j]], {i, 2, Length[partsS] - 1}, {j, 1,
       Length[partsS[[1]]]}]
     }^2]];

Fortunately, constrained optimization is built-in in Mathematica:

allVariables = Flatten[{variablesX, variablesY}];
constraints = -blockSize/3. < # < blockSize/3. & /@ allVariables;
initialValues = {#, 0} & /@ allVariables;
solution =
  FindMinimum[{matchEnergyFit + 0.1 smoothnessEnergy, constraints},
   initialValues];

Let's look at the result:

grid = Table[{(j - 0.5)*blockSize - dx[i, j], (i - 0.5)*blockSize -
      dy[i, j]}, {i, Length[partsS]}, {j, Length[partsS[[1]]]}] /.
   solution[[2]];
Show[src, Graphics[
  {Red,
   Line /@ grid,
   Line /@ Transpose[grid]
   }]]

The 0.1 factor before smoothnessEnergy is the relative weight the smoothness energy gets in relation to the image match energy term. These are results for different weights:

Possible improvements:

Like I said, perform a rigid registration first. With a white background, simple image moments-based registration should work fine.
This is only one step. You can use the offsets you found in one step and improve them in a second step, maybe with a smaller search window or smaller block sizes
I've read articles where they do this without blocks at all, but optimize an offset per pixel.
Try different smoothness functions


Determination of periodicity in data and finding mean I have to find whether there is any pattern (I mean periodicity or close to periodicity) and if there is, for one cycle i have to perform numerical integration to determine mean. In the first picture I have done moving average(every 254 data points) using excel.

Then in the next picture i have done again moving average(with 250 points) using the data got after performing first moving average.

I have also thought about doing Fourier transform.
These datas are values of cl (lift co-efficient) vs time in a unsteady solution of CFD simulation.Variation of values at initial time steps are not significant as some time is needed for simulation to give exact value.  DATA file is here (drive.google.com/open?id=0B1KEir1gFkaMQW5GbmFoZlpRUVk)

To determine periodicity, do an autocorrelation. The autocorrelation will have peaks at time offsets consistent with the periodicity.

You mentioned using the Fourier Transform, if convenient you can do a circular autocorrelation using the FT as follows:

$$Corr = ifft(fft(x) fft(x)^*)$$

Where $fft(x)^*$ is the complex conjugate (if you signal is complex).


Choices of convention and notation for the Fourier transform? The definitions of the Fourier transform and inverse Fourier transform I learned in college were

$$
F(j\omega) = \int_{-\infty}^{\infty} f(t) e^{-j\omega t}\ dt
$$
$$
f(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}F(j\omega)e^{j\omega t} d\omega
$$

The salient characteristics of this convention are

Non-unitary transform; frequency-domain units are radians (variable is $\omega$)
"Time-domain" units are in time (variable is $t$)
Function transforms are denoted by capital letters ($F$ vs. $f$)
The $j$ in $F(j\omega)$ strictly denotes that the function is a Fourier transform
And of course, the usual EE convention that $j=\sqrt{-1}$.

Nowadays I use a much different convention, essentially that used on
the wikipedias:

$$
\hat{f}(\xi) = \int_{-\infty}^{\infty} f(x) e^{-j2\pi\xi x}dx
$$
$$
f(x) = \int_{-\infty}^{\infty} \hat{f}(\xi) e^{j2\pi\xi x} d\xi
$$
The characteristics of this convention are

Unitary transform; frequency-domain units are normalized frequency (variable is $\xi$)
"Time-domain" units are unitless (variable is $x$)
Function transforms wear hats ($\hat{f}$ vs. $f$)
Variables in the Greek alphabet denote transformed variables in the Latin alphabet ($\xi$ vs. $x$)

I greatly prefer this convention for several reasons.

Using a unitary convention greatly increases the symmetry and clarity of Fourier duals:
compare

$\mathrm{rect}(x)\leftrightarrow\mathrm{sinc}(\xi)$,
$\mathrm{sinc}(x)\leftrightarrow\mathrm{rect}(\xi)$ to
$\mathrm{rect}(t)\leftrightarrow\mathrm{sinc}(\frac{\omega}{2\pi})$,
$\mathrm{sinc}(t)\leftrightarrow\mathrm{rect}(\frac{\omega}{2\pi})$.

Using $x$ instead of $t$ for the "time-domain" variable makes the equations more
agnostic of one's problem domain. This makes it much easier to analogize about
2D image processing concepts in terms of 1D signal processing concepts, without
the cognitive dissonance of using $t$ as a variable representing
distance, or having to change variables around when going from one domain to the other.
I find capital letters to be more useful for denoting discrete-valued variables/functions than for representing transformed functions.
Using the hat more clearly denotes the Fourier transform as an operator that is
applied to $f$, where the resulting function accepts a frequency-domain parameter $\xi$. Compare this to the far more ungainly $\mathcal{F}\{f\}$, the
"traditional" way I learned in college to denote the Fourier transform as an operator,
which just seemed way too confusing to me at the time ($\mathcal{F}{f(t)}$ vs
$\mathcal{F}\{f\}(t)$ vs $\mathcal{F}\{f\}(\omega)$ etc.)
I generally found that doing signal processing analysis with radians just sprinkled
a lot more $\pi$s around than I felt necessary. Using units of normalized frequency
makes a lot more sense to me, particularly when working through problems involving
sampling theory.

Of course, it would be quite vain of me to consider my choice of convention to be
superior to that used by others. But I'm having a hard time coming up with good
reasons to prefer the convention I originally learned in college (ie, reasons
that don't involve tradition).

Presently, I can think of one decent reasons for preferring the "traditional" convention:
Using the non-unitary transform, and the $F(j\omega)$ parameter notation, greatly
improves notaional consistency with the
Laplace transform. Also, hats might be easier to lose/confuse than capital letters.

Can anybody think of other reasons to prefer the "traditional" (non-unitary) convention?
Is this "traditional" convention the same as what you learned a signal processing course
(if you took one)?
Which convention do you prefer?

The choice of convention should be the one most appropriate (or familiar) for the audience you are trying to communicate with.


Best (Perceptually / Objectively) Super Resolution Methods Out There? I'm curious about the advances in the area of image super-resolution (SR) that have given the best results to date, both perceptually (visually pleasing) and objectively (e.g. PSNR, SSIM criterias). I'm mostly talking about Single Image SR (SISR), but I'm open to input on Multiple Image SR as well.

From what I've researched, Photo-Realistic Single Image Super-Resolution Using a Generative Adversarial Network looks like the one, at least as far as SISR goes. However, I am unsure as to what alternatives are out there, that come close or maybe even surpass it (both in the field of neural networks and other, more traditional fields).

Remark
If this reaches anyone's interest but my own, perhaps this can be a question that gets new answers as time and advances move forward.

In our days the Deep Neural Network methods certainly are generating best results.
Due to the intense research going on in this field the best method is a moving target hence one can not pin point to one.
One generation before them the best methods were based on Dictionary Learning.
For example you can use the K-SVD for Single Image Super Resolution.
Those approaches / methods works pretty well and are even improved using the EPLL approach.
You should be aware those methods are really slow for real world images.
Probably you should invest in the DNN world.
Pick any one of the latest on ArXiv Sanity.


how could correlation-based image patches matching works? In general, correlation-based image matching can be referred to as
$$c(I_1,I_2) = \sum_{x \in patch} f(I_1(x), I_2(x))$$
, I do think normalized cross-correlation like
$$ncc(I_1, I_2) = \sum_{x \in patch}\frac{(I_1(x) - \mu_1)(I_2(x) - \mu_2)}{\delta_1\delta_2}$$.
is a good indicator to show if 2 patches match or not.

BUT my problem is that, if
$$\begin{eqnarray*}
f(I_1(x), I_2(x)) &=& \sum_{x \in patch}I_1(x) \cdot I_2(x)
\end{eqnarray*}
$$
, then I don't think $c(I_1, I_2)$ will be a good, right?

Furthermore, this type of correlation-based matching techniques won't work well if one of the images rotates, right?

Normalized is preferred to regular correlation to handle variation in brightness, or intensity, between images if nothing else.

Regarding your second comment, the correlation is going to be lower if the shape rotates (unless it happens to be symmetrical in all four directions).  There are other features you can consider that are rotation invariant as well.  See this post for some examples.


In speech recognition applications, how are different voice amplitudes handled? (i.e loud vs. quiet tones) I am doing some research related to real-time voice/hotword recognition engines. In most current implementations, input is divided into frames (overlapping or not), and audio features are extracted per frame (most common being MFCCs) and fed into a Hidden Markov Model or a Neural Network of sorts.

Most papers I read address issues such as noise removal/reduction (using methods like Cepstral Mean Normalization), however I couldn't find any mention of how different voice amplitudes are handled. For example, I can train an engine to recognize my voice when I speak normally, however if I change my volume (speak louder), then the extracted features would look different (same shape but larger magnitude). Since this is a real-time system, I am not sure how real-time normalization is even possible, or whether it should be applied on the voice samples, or the extracted features. Or perhaps it is solved by training the system on variations in the speaker's volume?

Your help is greatly appreciated.

Some of the channel effects can be indeed removed by doing the Cepstral Mean Subtraction/Normalization. Nevertheless that generally applies only to "convolutive" distortions that are constant. Any additive distortions, i.e. white noise, babble noise usually cannot be removed via CMS. But like you said this topic is handled via different methods to cope with the noise removal.

As for the immunity of MFCC's with respect to different audio levels, you can do it easily by taking the appropriate coefficients. Generally the first MFCC coefficients is obtained by fitting the constant value curve ($\cos(0)$) to your log-energy filter banks. Therefore it is highly correlated to the RMS energy of your signal. If your remove that coefficient (often called 'static') then in theory you make your model volume (gain) independent. The rest of coefficients is not really related to the energy of your signal. Usually researchers drop the first coefficient, but they are adding it first and second derivative to the $\Delta$ and $\Delta\Delta$ features.

Obviously there is more to that. For example the Lombard Effect might take place and it's obviously changing the envelope of your signal.


Reference for Supersampling I want to downsample images to arbitrary sizes using supersampling to avoid aliasing effect.

The only two good explanations I found were on Wikipedia and everything2.com, but there are still gaps. For example:

The image should be upsampled before supersampling. But how much? Is there any special filter to be used for upsampling prior to supersampling? Is just taking interpolated color values on sub-pixel positions in source image sufficient?
What if I want to resize image to e.g. 80% of its size with supersampling? How much I have to upsample it and how to treat non-integral size ratio (e.g. downsampling from 379 pixels to 377 pixels)?
Are the multiple samples taken from the single pixel, multiple pixels or around sample point (both cases can occur)?

UPDATE:

I have tested the "Super Sampling" method in Paint.NET with test target on this page. Suprisingly, the result looked just like "Photoshop Bicubic" filter:

The only advantage can been achieved by pre-blurring input image and than using some conventional method.

I've implemented and tried ordinary "Box" and "Cubic" convolution filters with excellent results!

So now I don't see any benefit of using supersampling over the convolution-based filters. Or is there any?

EDIT: I took this answer down for a time because I realized a lot depends on whether there are any idiosyncrasies in how images are mapped to pixels by actual devices, and, not being an image guy, I don't know a lot about that.  I decided to bring it back, with the caveat that the answer may be insufficient given said idiosyncrasies.

  The image should be upsampled before supersampling. But how much?

You are just trying to soften quantization (pixel resolution) error by collecting information from neighboring pixels.  Thus, oversampling by a factor of 2 should be plenty.  In other words, once you oversample until your sample resolution is twice that of the pixel resolution, nothing can be gained by oversampling more.

  Is there any special filter to be used for upsampling prior to supersampling?

I might be misunderstanding the terms, but it appears to me that upsampling is part of supersampling.  No, you don't need to use any special filters.  Insert 0's in between your samples and then low-pass filter.  You'll probably need to adjust the gain of your low-pass filter by your upsampling rate to get the pixel intensity to come out right.

  What if I want to resize image to e.g. 80% of its size with supersampling?

Then factor that into your upsample/decimation calculations.  After you upsample you have to decimate down to the pixel resolution.  By resizing the image you are, in effect, changing the pixel resolution.

  How much I have to upsample it and how to treat non-integral size ratio?

How much you have to upsample it depends on the details of the image, pixel resolution, etc.

  How much I have to upsample it and how to treat non-integral size ratio?

You can deal with non-integral ratios with fractional resampling.  Essentially you upsample by an integer and decimate by another integer, giving you a fractional sample change.


Radon Transform reducing the number of parallel beams (MATLAB) I am trying to understand different filters' and other things' effects on Radon Transform by using MATLAB. First I upload an image, then take its Radon transform with radon function and then reconstruct it with iradon function.

Now I wanted to see the difference between undersampling the projection data by reducing the number of parallel beams by a factor of 4 and reducing the number of rotation angles by a factor of 4. However according to what I understood from radon and iradon functions' descriptions, I can only change the number of rotation angles. Is not it true? If not, what is the way to do it?

Here is what I do to decrease the number of rotation angles by a factor of 4:

I = imread ('asdasd.tif');
Rad = radon (img,0:4:179);
I_ram = iradon(Rad,0:4:179,'linear','Ram-Lak');
I_ham = iradon(Rad,0:4:179,'linear','Hamming');

Unfortunately you can't do this with the radon function as it comes.  However, the function is written in MATLAB form (as opposed to being a compiled function), so you could extend the function yourself.  This would also be a good way to learn more about the process.

To reconstruct properly, you would also need to perform edits on the iradon function.

An alternative way to get what you're looking for would be to get projections as they come from radon, then reduce the array size by averaging the projection rays in groups of 4.  You would probably want the reconstruction to be the same size, so you'd need to repeat the averaged values 4 times in the argument to iradon.  So maybe something like this:

p = phantom ();
proj = radon (p, 0:179);
proj_avg = proj (1:4:end,:) + proj (2:4:end,:) ...
         + proj (3:4:end,:) + proj (4:4:end,:);
proj_rep = reshape (repmat (proj_avg',1,4)', [], length (proj(:,1)))';
recon_few_beams = iradon (proj_rep, 0:179);


Do antisymmetric FIR filters preserve the shape of their input waveforms? For antisymmetric FIR windows, the phase response is of the form $a\omega+b$, where $b$ is $\pm\frac{\pi}{2}$, so it is not strictly linear. Does the waveform of the signal gets retained?

Such phase response is called Generalized linear phase, where you are allowed the $\pm \pi/2$ constant phase, and zero crossings (which add $\pi$ jumps in phase at one or more frequencies).

Filters with an antisymmetric impulse response do in fact have the $\pm \pi/2$ constant phase.

Simplest case: $h[n] = \delta[n] - \delta[n-1]$, with frequency response $H(e^{j\theta}) = 2j e^{-j \theta /2} \sin( \theta/2 )$.

In practical terms, these filters are types of differentiators: they all have a zero at $z=1$ (0 frequency) due to the fact that $\displaystyle\sum h[n] = 0$.

It is clear that these filters do not retain the waveform of the signal.

Normally you say that linear phase filters "retain the waveform" in case of low-pass filters, where you want to filter out noise but retain the shape of a pulse. More importantly, you want a good eye-pattern when decoding, and a non-linear-phase filter will result in asymmetric step responses and a poor eye pattern.


Simple fft to Gaussian pulse with MATLAB trying simply to create femtosecond pulse in MATLAB exactly like in the attached image,

the carrier frequency is around 374THz, and my sampling frequency is 10 times the carrier.
The results of the fft yields nothing understandable... tried to change some of the variables, fs,t-around zero, fftshift...
but could not conclude what's wrong with code.

my code follows matlab fft example :

clear all ; close all ; clc
f=374.7e12;%Thz
fs=f*10; %sampling frequency
T=1/fs;
L=1000;
sigma=5e-15;
t=(0:L-1)*T; %time base

x=(exp(-(t-50e-15).^2/(2*sigma)^2)).*exp(-1i*2*pi*f*t);
subplot(2,1,1)
plot(t,real(x),'b');
title(['Gaussian Pulse \sigma=', num2str(sigma),'s']);
xlabel('Time(s)');
ylabel('Amplitude');
ylim([-1 1])
xlim([10e-15 90e-15])

NFFT = 2^nextpow2(L);
X = fft(x,NFFT)/L;
Pxx=X.*conj(X)/(NFFT*NFFT); %computing power with proper scaling
f = fs/2*linspace(0,1,NFFT/2+1); %Frequency Vector

subplot(2,1,2)
plot(f,2*abs(X(1:NFFT/2+1)))
title('Magnitude of FFT');
xlabel('Frequency (Hz)')
ylabel('Magnitude |X(f)|');

the results are :

You simply don't plot what you want to see. Note that your time domain signal is complex-valued and you modulate by a negative frequency. So the range of frequencies where things are happening are the negative frequencies (or - by periodicity - the frequencies in the range $[f_s/2,f_s]$), but those you don't plot. If you change the definition of your frequency vector and the corresponding plot command you'll see what you expect to see (i.e. a Gaussian in the frequency domain centered at $f_s-f=3.37e15$ Hz):

f = fs*linspace(0,1,NFFT); % (full range) Frequency Vector
...
plot(f,2*abs(X))

OR, simply change the negative modulation frequency to a positive one

x=(exp(-(t-50e-15).^2/(2*sigma)^2)).*exp(1i*2*pi*f*t);

and leave everything else unchanged. Then the Gaussian in the frequency domain is centered at $f=374.7$ THz. This is probably what you expected to see in the first place.


What type of correlation is this equation? In a book [1], gives an equation (equation 5.31 in the book) for estimating frequency offset in OFDM systems proposed by Classen. the procedure is:

  first, two OFDM symbols, $y_l [n]$ and $y_{l+D} [n]$ , are saved in a memory after synchronization. then, the signals are transformed into ${Y_l [k]}_{k=0}^{N-1}$ and ${Y_{l+D} [k]}_{k=0}^{N-1}$ via FFT, from which pilot tones are extracted. After estimating CFO from pilot tones in the frequency domain, the signal is compensated in the time domain.

Here is the equation:

$ \hat \epsilon_{acq} = \dfrac{1}{2\pi \cdot T_{sub}} \underset{\epsilon}{max} \left \{ \left \vert \sum_{j=0}^{L-1} Y_{l+D}[p[j],\epsilon]\cdot Y_{l}^{*}[p[j],\epsilon]\cdot X_{l+D}^{*}[p[j]]\cdot X_l[p[j]] \right \vert \right \}  $

where $L$, $p[j]$, and $X_{l}[p[j]]$ are the number of pilot tones, the location of the jth pilot tone, and the pilot tone located at $p[j]$ in the frequency domain at the $l^{th}$ symbol period, respectively.

From the equation, It's supposed to find the normalized CFO, $\epsilon$. but, the parameter we are finding is involved in the equation $[P[j],ε]$.
My question is: how this parameter we're supposed to estimate is in the middle of the equation? What type of correlation is this? Or did I misunderstand the equation? Please explain.

[1] Yong Soo Cho et al., "MIMO-OFDM wireless communications with MATLAB", John Wiley & Sons (asia), 2010

I understand your confusion because the equation is barely understandable from the information given in the book. It becomes more clear from the original paper by Classen and Meyr [1] from which it has been taken. They propose a two stage frequency offset estimation that consists of an acquistion stage and a tracking stage. The equation you've cited represents the acquisition algorithm. In fact, it is a searching algorithm that "tries out" various values of frequency offset $\epsilon$ and chooses the one value $\hat\epsilon_\mathrm{acq}$ that yields the maximum correlation sum. It can be described as follows

For every $\epsilon$

Compensate the received time domain signal for the supposed frequency offset $\epsilon$
Calculate the DFT of the frequency offset compensated RX signal. This yields $Y_{l}[p[j],\epsilon]$
Calculate the correlation sum $R_\epsilon$. (The term inside $|\cdots|$)

Then, find the maximum $\hat R_\epsilon$ of all $R_\epsilon$ and calculate the estimated frequency offset by
$$
\hat \epsilon_{acq} = \dfrac{\hat R_\epsilon}{2\pi \cdot T_{sub}}
$$

Of course it's not feasible to evaluate the correlation sum for every $\epsilon$ because the number of possbile values is infinite. Instead, you have to restrict the range of values to a reasonable frequency offset that is determined by the hardware you're using. Additionaly, this "search range" has to be partitioned into discrete values. The authors of [1] note

  In practice we found that it is sufficient to space the trial parameters $0.1/T_{sub}$ apart from each other.

Where $T_\mathrm{sub}$ is the length of one OFDM symbol including guard interval.

[1] Classen, F. and Myer, H. (June 1994) Frequency synchronization algorithm for OFDM systems suitable for
communication over frequency selective fading channels. IEEE VTC’94, pp. 1655 1659.


Errata for Lyons's DSP book After reading the comments$^1$ and looking through various texts, I am going to use Understanding Digital Signal Processing for my self-study (Oppenheim & Wilsky was my 2nd choice). I'm 99% sure that a solutions manual doesn't exist. But does anyone know where I can find the errata for the 3rd edition?

I saw an old thread where the author, Dr. Richard Lyons offered to send it to someone so I can email him personally, assuming I can find his address. Before going that route, I figured that someone here may have it. Thanks.

$^1$${\tiny\textrm{ Thanks to all for past threads where this was discussed and great advice given.}}$

This blog post by Richard Lyons contains the links to the Errata for all editions of Understanding DSP

International 3rd Edition Errata here.


Power spectral density vs. FFT bin magnitude What's the difference between these? Both are measurements of some form of signal power, but surely there's some difference between the power they are measuring?

The fast Fourier transform ($\textrm{FFT}$) algorithms are fast algorithms for computing the discrete Fourier transform ($\textrm{DFT}$). This is achieved by successive decomposition of the $N$-point $\textrm{DFT}$ into smaller-block $\textrm{DFT}$, and taking advantage of periodicity and symmetry.

Now, the $N$-point $\textrm{DFT}$ of a sequence $\{x[0], x[1],\cdots, x[N-1]\}$ is:
\begin{equation}
X(f_k) = \displaystyle \sum_{n = 0}^{N - 1} x[n]\exp\left(-j2\pi f_kn\right) \tag{1}
\end{equation}
For $f_k = k/N$ and $k = 0, 1, \cdots, N - 1$. And the $\textrm{FFT}$ magnitude at bin $k$ is the $\textrm{DFT}$magnitude at bin $k$. For a given $N$ that is:

\begin{equation}
\left|X(f_k)\right| = \left|X\left(\frac{k}{N}\right)\right| = \left|X(k)\right| = \displaystyle \left|\sum_{n = 0}^{N - 1} x[n]\exp\left(-j2\pi nk/N\right)\right| \tag{2}
\end{equation}

You talk about power spectral estimation when the signals being analyzed are characterized as random processes. With random fluctuations in such signals, statistical characteristics and average characteristics are normally adopted. For a wide sense stationary $(\textrm{WSS})$ discrete random process, the PSD is defined as:

\begin{equation}
P(f) = \displaystyle \sum_{m = -\infty}^{\infty} r_{xx}[m]\exp\left(-j2\pi fm\right) \tag{3}
\end{equation}

For reasons you can find in this answer, you see that the squared magnitude of the signal's $\textrm{DFT}$ is taken as the estimate of the PSD in most practical situations. One form, among other variations/methods, is:

\begin{equation}
P(f_k) = \frac{1}{N} \displaystyle \left| \sum_{n = 0}^{N - 1} x[n]\exp\left(-j2\pi f_kn\right) \right|^2 \tag{4}
\end{equation}

  What's the difference between these?

Comparing $(2)$ and $(4)$, you have:

\begin{equation}
P(f_k) = \frac{1}{N} \left| X(f_k)\right|^2
\end{equation}

From the bin number $k$ to frequency in $\textrm{Hz}$, $F = \frac{F_s}{N}k$

For more reading on PSD estimation check this question, that question, and this question.

EDIT:

The power spectral density, $\textrm{PSD}$, describes how the power of your signal is distributed over frequency whilst the $\textrm{DFT}$ shows the spectral content of your signal, the amplitude and phase of harmonics in your signal. You pick one or the other depending on what you want to observe/analyze. And no they're not the same as you can see from the equations above and links given. Their spectra are generally not the same. One is estimated as the squared magnitude of the other.


How to remove sum and frequency components in FFT spectrum during signal resampling.? I am working on to remove non linear tuning of laser . I am using resampling technique based on triggering technique .
Firstly, I make collect a vector of maxima and minima from Tigger sine wave of frequency F1. Then I reconstruct signal sine wave (of frequency F2) with the maxima and minima vector from the trigger. FFT is performed after this step to identify noise free frequency components of the signal.
This kind of of resampling is working excellent. But after FFT,  besides the actual frequency component of the signal wave, it also gives rise other unwanted frequency components.
These are exactly sum and frequency components arising from F1 and F1,( i.e., F1+F2 & F1-F2 ).
For the clear understanding of my question please refer to attached image  and suggest me the possible solution to avoid these unwanted frequency generations.
thank you

You should be using an appropriate low pass filter during reconstruction.


Finding Accurate BPM and Beat Time Values in Audio I working on writing a program to get the BPM rate and the Beat Times in audio songs but I am having trouble trying to come up with reliable testing strategies for my design.What do I benchmark my BPM rates and Beat Times against? Should I just use other softwares or is there ways of checking if my BPM rates and beat times are correct? (apart from physically listening and tapping along to the song).

Comparing your results against data annotated by a human listener is the way to go. You can try downloading some of the evaluation data from MIREX or look for other research databases.


Attempting to reverse mid conversion I'm new to to DSP in general and have been trying my hand at making some very simple VST's just for practice purposes.

I am trying to convert the stereo channels into Mid (in the left channel) and Side (in the right) without using a placeholder to store output. Here is a very basic version of the code I am currently trying:

    *LeftSample = (*LeftSample + *RightSample) * 0.5;
    *RightSample = (((*LeftSample * 2) - *RightSample) - *RightSample) * 0.5;

The first line assigns the mid output to the left channel, but in the second line I am trying to reverse the mid conversion so that I can use the left channel's original output when performing the L-R operation to create the side channel in the right output.

Am I missing something in thinking that this should be achievable? Is there some reason that

    (*LeftSample * 2) - *RightSample

does not reverse the mid conversion, (L+R) * 0.5? Is it not possible to derive the original left channel using the mid signal and the right signal?

Again, I'm new to this, so I apologize for any errors in terminology I have made, and for anything that has been poorly communicated. I appreciate all feedback and criticism (please, tear any errors I've made to pieces). I've looked around quite a bit, and haven't had any luck; if anyone has any resources or directions to point me in, I would greatly, greatly appreciate it!

You can fully reconstruct L and R from M and S, as you suspect. In your case:

$ M = {{L + R} \over 2} $

$ L = 2M - R $

$ S = {{L - R} \over 2} $

then

$ S = {{2M - R - R} \over 2} $.

So the math checks out. Are you not getting the results you expect?

When I've worked with M&S in the past, I've alyways just done the scaling by 1/2 when returning to L&R. I guess it depends on what you're intending on doing with the M&S signals though. If you want the sum of M and S to have equal energy to the sum of L and R, you should scale by $ 1 \over \sqrt 2 $ when converting both to and from M&S.


Plotting a histogram of particle intensity vs pixel I am trying to reproduce a figure, shown below, of a particle intensity binned by the pixel number then fit the distribution to a Gaussian. I have a similarly magnified intensity image of my particle, with 8-bit depth, loaded into MATLAB.

I am having difficulty understanding how the bottom figure seems to retain spatial information after binning. The goal is to fit the distribution to a Gaussian.

The expected distribution of intensity versus x,y position is formed by two idependent and identically distributed gaussian random variables; one for each axis. The distribution shown is for one "slice" through this 2D distribution.  The top view of the distribution can be envisioned as a 3d Gaussian "hill", and the histogram shown is for one slice through this hill as viewed from the side, where the slice can be from any direction as long as it goes through the center (origin).


Compute filter output with given values of $n$ \begin{align}
x[n] &=\begin{cases}
 1& \text{for}\quad n \ge 0\\
0& \text{otherwise}
\end{cases}\\
h[n] &= \delta[n] - \delta[n-1]\\
y[n] &= x[n]*h[n]
\end{align}

  Find the values of $y[-1]$, $y[0]$, $y[1]$, $y[2]$.

I thought $\delta[0] = 1$ and $\delta[n] = 0$ otherwise, and that $y[n] = x[k]h[n-k]$:

\begin{array}{|c|c|c|c|c|c|}\hline
n&x[n]&d[n]&d[n-1]&h[n]&y[n]\\\hline
-1& 0 & 0 & 0 & 0&0\\\hline
0&1&1&0&1&1\\\hline
1& 1 & 0 & 1 & -1&-1\\\hline
2& 1 & 0 & 0 & 0&0\\\hline
\end{array}

So I put $\text{ 0  1  -1 0 }$ as my answer but apparently that's wrong. What did I mess up? Did I misinterpret the meaning of $\delta[n]$?

Your impulse response is $h[n]=\delta [n] - \delta [n-1]$, so it basically takes a copy of the input and then substracts to it a delayed version of the former. When $x[n]=u[n]$, the output signal is:

$$y[n] = u[n] - u[n-1]$$

This means that we have a step function (i.e. it is equal to $1$ for $n \geq 0$ and $0$ otherwise), and we substract to it a shifted step (i.e. it is equal to $1$ for $n \geq 1$ and $0$ otherwise). Therefore, $\forall n \geq 1$ we are cancelling out the two signals ($1-1=0$). This leaves us with only one not-zero value, which happens to be at the origin. So the output turns out to be

$$y[n]=\left\{
 \begin{array}{ll}
  1  & \mbox{if } n = 0 \\
  0 & \mbox{otherwise }
 \end{array}
\right. \implies y[n]=\delta[n]$$


How to change speed of audio samples without changing pitch? This question may be more appropriate on music stack exchange but since sound is a branch of physics and this is more of a technical question relating to sound waves than a musical question I asked it here.

How can you change the speed of audio samples without changing pitch?   I have heard songs that are sped up without changing the pitch of the song.   I also know that simply compressing an audio sample changes the frequency and thus the pitch of the sound waves that are perceived by our brains.   So how is this done without changing pitch?

This question (about "time-scaling" audio) is closely related to pitch shifting, which is time-scaling combined with resampling.  But changing the speed without changing pitch is only time-scaling, so there is no resampling involved (contrary to what thomas has suggested).

There are frequency-domain methods (phase-vocoder and sinusoidal modeling) that can change speed for an orchestral mix (or some other broadbanded non-periodic sound) without glitches, but they can become a little "phasey" in their stretched sound.

If it's monophonic, like a solo or a single note or string or voice, a less expensive time-domain approach is sufficient and can sound very good, under this monophonic condition.  If it's a single note or tone, the signal is quasiperiodic.  That means, in the vicinity of some time, $t_0$, that

$$ x(t) \approx x(t - \tau(t_0)) \qquad \text{for } t \approx t_0 $$

So $\tau(t_0) > 0$ is the estimated period of the quasiperiodic waveform in the neighborhood of $t_0$.  This estimated period $\tau(t_0)$ is estimated by finding the value of $\tau$ that minimizes something like

$$ Q_x(\tau, t_0) = \int \Big|x(t) - x(t-\tau)\Big|^2 w(t-t_0) \, dt $$

where the window $w(t-t_0) \ge 0$ is centered around the time $t_0$.  If there are multiple candidates of $\tau$ that minimize $Q_x(\tau, t_0)$, usually it is best to pick the smallest $\tau$ in which $Q_x(\tau, t_0)$ is small.

Now imagine copying your note $x(t)$ and offsetting (delaying) the copy in time by the amount $\tau(t_0)$.  That delayed copy would be $x(t - \tau(t_0))$.  Then you have two identical waveforms, except for the offset of $\tau(t_0)$ but around the time $t_0$ the two waveforms will look almost the same because the offset is exactly one period, based on the estimate of the period around time $t_0$.

Then you can cross-fade (using a Smoothstep function $S(t)$) from the original to the offset copy, the cross-fade will not have any nasty cancellation or "destructive interference" because the waveforms are lined up.  And the result will be the same note, but one period longer.  (Longer by $\tau(t_0)$ seconds.)

$$ y(t) = \Big(1-S\big(\tfrac{t-t_0}{\tau(t_o)}\big) \Big) \, x(t) \ + \ S\big(\tfrac{t-t_0}{\tau(t_o)}\big) \, x(t - \tau(t_0)) $$

(That won't stretch it long enough to notice a difference, but if you repeat this operation many times per second, you can stretch a patch of sound by even a factor of two.)


Unstable plant transfer function identification I would like to find the transfer function of an unknown unstable SISO plant.

If it was a stable plant, I would input a sine sweep and measure the frequency response at the output; but I cannot do this since the plant is unstable.

I use a controller to get a stable closed-loop.

Methods I have tried:

Sine-sweep at the input of the stable closed-loop. Problem: I have trouble extracting only the plant (without the controller) from the closed-loop transfer function.
Add a sine disturbance (at different frequencies) at the controller output, and then estimate the gain of the transfer function from the input of the plant to the measurements of the plant by dividing the peak-to-peak of each signal, and the phase by calculating the time between the input's peak to the measurements's peak (for each frequency).

This method takes a lot of time since it requires re-tuning the controller gains for different frequencies. Otherwise, with high gains at low frequencies, the input disturbance is attenuated, or with low gains at high frequencies, the output is saturated/unstable.

Also, this method gives periodic but non-sinusoidal signals. I'm not sure I'm calculating the gain and phase of the plant transfer function correctly in this case.

I have two questions:

Is the second method correct?
Is there a better way to find the transfer function of an unstable plant?

See the image (found in this article).

Inject a signal into the loop somewhere (I show it below as going between the controller and the plant).  Then measure the input to the plant and the output from the plant.  Note that the image is a bit confusing -- it assumes that the normal command to the loop is 0 or some constant.  It would be better to actually measure right at the output of the plant.

Then use the input to and output from the plant to estimate the plant response.

You can do this for swept sine, impulses, random noise, and steps.  If you're concerned about noise, you can take long samples -- when I do swept-sine measurements using this technique, I set things up so that I can take arbitrarily long measurements at each frequency, to make sure I can capture the real signal in the presence of noise.  If you're doing steps, you can get the response to an arbitrarily high number of steps, average, and then do ARMA or whatever method you wish to extract the plant transfer function.  Ditto random data or impulses.


The z factor in polyphase decomposition I'm a beginner in DSP and I'm going through the textbook of Oppenheim's Discrete Time Signal Processing. There are two figures in the text, one which I can visualize, and the other I can't. The first figure shows the polyphase decomposition of an impulse response function:

I get why the $z^{-n}$ is placed after the upsampling operation. After the zeros are added during the upsampling operation, appropriate delays need to be added so that the sum gives the original impulse response. The next step taken in the textbook isn't clear to me. I do not get why the $z$ delays are added before the polyphase components in the second figure, as shown below:

Doing so suggests the delay operation is commutative with the expansion, but I cannot convince myself that it is right. I tried doing that for the first figure and I do not get the same results. Is there something I'm misunderstanding here?

Thanks.

The polyphase components of $h[n]$ are

$$e_k[n]=h[Mn+k]\tag{1}$$

Note that if $h[n]$ starts at $n=0$, then so do all polyphase components $e_k[n]$. In order to get back the original sequence $h[n]$ from the polyphase components $e_k[n]$, we need to insert $M-1$ zeros between the samples of $e_k[n]$, and shift each component back to its original position, i.e. shift $e_k[n]$ by $k$ samples. In the frequency domain this corresponds to

$$H(z)=\sum_{k=0}^{M-1}E_k(z^M)z^{-k}\tag{2}$$

Filtering an input signal with $\mathcal{Z}$-transform $X(z)$ results in

$$Y(z)=H(z)X(z)=\sum_{k=0}^{M-1}X(z)E_k(z^M)z^{-k}\tag{3}$$

which corresponds to the structure in figure $4.37$ in Oppenheim and Schafer. The terms $E_k(z^M)z^{-k}$ can clearly be implemented either by delaying the input or the output of the polyphase components:

$$X(z)\left[E_k(z^M)z^{-k}\right]=\left[X(z)z^{-k}\right]E_k(z^M)$$


Wavelet transform and FFT using to extract feature power bands with EEG signals I am using 5 channels [ fz , cz , c3 , c4 , pz] to detect drowsiness of driver
My First Question is, what is the right input to get feature power band ( Theta , alpha , gamma , beta ) to wavelet transform ? ( these 5 channels or 1 channel or what ? )
My Second Question is, Is it right to classify data based on theta only got from wavelet transform ?

Alpha rhythm
Generally, alpha-band oscillatory activations (8-10Hz) relate to relaxation and in principle accompanied with closure of eyes. This is the prime marker that is used to detect drowsiness, but surely not the only one (see alpha dropout, NREM1, eye-rolling EEG artefacts).
Event detection
To detect alpha, electrodes from occipital or adjacent regions are used. Due to low SNR in EEG signals and its variance across EEG systems, you may want to try different detection methods to decide for your classifier.


Using FFT as a channeliser This article mentions that the DFT/FFT can be used as a channeliser in a similar way to a filterbank. I get how the DFT filters into spectral bins, but its output is in the frequency domain, i.e. Fourier coefficients, nolonger in the time domain, so either I'm missing something or they've neglected to mention that DFT filterbank output would need to be transformed back to the time domain to function like an actual filterbank. Can someone shed light on this?

First, I suggest you take a look at the STFT.

Let say you use $N$-point FFT, then you filterbank has $N/2$ filters (assuming real signal). Usually, you divide your signal into $K$ segments, each $N$ samples long (could be overlapping) and apply the FFT to each segment. Every time you apply the FFT you get a single sample for each filter. Therefore your filterbank has $N/2$ filters with $K$ samples.


Effect in Fourier-domain when the signal is "zero-stuffed" in time-domain Suppose $x[n]=[x_0,x_1, \ldots,x_{N-1}]$ and $y[n]=[x_0,0,x_1,0,\ldots,0,x_{N-1}]$. What is the relationship between the Fourier Transforms of both the signals?

When I try for some specific examples in MATLAB, I see that $Y[k]$ has an image of $X[k]$ but how do we show this mathematically?

Let's be more general and consider the insertion of $L-1$ zeroes between each sample. Mathematically, this can be written as
$$
y[m]=
\begin{cases}
x[n] & \text{if } m = nL\\
0  & \text{otherwise}
\end{cases}.
$$
The DTFT of $y[n]$ is then given by
$$Y(e^{j\Omega}) = \sum_{m=-\infty}^{\infty} y[m]e^{-jm\Omega}.$$
Because $y[m]$ is non-zero only when $m$ is a multiple of $L$, that is $m = nL$, this becomes
$$Y(e^{j\Omega}) = \sum_{n=-\infty}^{\infty} x[n]e^{-jnL\Omega} = X(e^{jL\Omega}).$$
From this, you can see that the DTFT of $y$ is the DTFT of $x$ with a period $L$ times smaller (i.e. $2\pi/L$ instead of $2\pi$). Intuitively, this makes sense: you know that "expanding" a signal in the time-domain is equivalent to "squeezing" its spectrum in the frequency domain (see the time scaling property of the Fourier transform).


Output length of an IIR Filter In general, the output length of an FIR filter is given by (N+M-1) samples, where N stands for signal length and M stands for the number of filter coefficients. Can the same formulae be applied for IIR Filters? Let's say I have 4 input samples, 3 feed-forward coefficients, and 2 feedback coefficients. Then can I tell that my output length is 4+(3+2)-1? It is also known that the impulse response of an IIR Filter is not of finite duration. Does this thing affect the output length of an IIR filter, in any way?

To answer your last question, the impulse response having infinite length affects the output length in every possible way. In general, the output sequence of an IIR filter has an infinite length, even for a finite length input sequence. Of course, for stable filters and finite-length input sequences, the magnitude of the output samples decreases, and the output values become virtually zero after a while.

There are, however, specific input sequences that result in cancellations in the output, and, consequently, the corresponding output sequences have only a finite number of non-zero elements. A simple example is the first-order system

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

You can verify yourself that for the specific input sequence $x[n]=\delta[n]-a\delta[n-1]$, the output sequence is finite and simply equals $y[n]=\delta[n]$. However, note that the use of finite precision arithmetic in an actual implementation may still result in non-zero output samples, even in cases when ideally the output samples should cancel.


Why is a negative exponent present in Fourier and Laplace transform? could anyone explain why there is a need of negative exponent in fourier and laplace transform.I looked through the web but I couldn't get anything.Does anything happen if a positive exponent is placed in these transforms.

While looking through http://1drv.ms/1tbV45S it says that if $s>0$  it becomes a rapidly decreasing function while if $s<0$ it becomes an rapidly increasing functin of t.I couldn't understand that.Can anyone illustrate this.

Matt is correct that the sign is convention.  I think that there is a reason for it beyond that though.

If we look at complex frequencies in the complex plane, they look like a constant vectors that rotate in one direction or another.  Positive frequencies rotate counter-clockwise, negative frequencies rotate clockwise, and "0 Hz" frequencies don't rotate at all.

The Fourier transform has a negative sign to intentionally rotate in the opposite direction as the frequencies that they are "looking" for.

The reason for the opposite rotation is that when the two frequency vectors are multiplied, their phases will repeatedly cancel out, so when the results are summed together there will be a massive vector due to all of the individual vectors lining up.

$$
X(f) = \sum\limits_{n=0}^{N-1}x(n)e^{-j2\pi kn/N}
$$

This is how the Fourier transform "looks" for frequencies.  If the two frequencies are the same or "close" (how close they need to be depends on the length of the DFT) they will line up well and cause a massive response in the summation.  I have showed how this works for the discrete Fourier transform (DFT), but the exact same reasoning applies to the continuous transform.

Hopefully this explains why the Fourier transform would want the vectors rotating in the opposite direction.  To be perfectly honest I don't know the Laplace transform well enough to give solid reasoning for its negative sign.  Since the two transforms are very closely related though (the Laplace transform being a generalization of the Fourier transform), I assume that it is for similar reasons.


Is it possible to find precise peak frequency using FFT phase shifts? Been experimenting with FFT on a generated sinusoid and found something strange that doesn't seem to be described anywhere (though I may be missing something of course).

A sinusoid that exactly corresponds to a bin gives what we all know: an amplitude and a phase shift of the sinusoid in that bin.

However, for a frequency that falls between two bins the phase shift of the strongest bin expressed in polar angle seems to be shifted by almost exactly the amount by which the peak frequency is shifted from the strongest bin number. The phase shifts of the neighbouring bins also seem to show some correlation but it is not as strong.

For example, if the distance between bins is $\Delta = 44Hz$, I generate $f = 440Hz$, I get a clear sinusoid at bin $i = 10$ that yields phase angle $\phi = 0$ (the sinusoid itself is shifted by a quarter period for simplicity, hence $0$).

If I generate say

$f = 440Hz + k·\Delta$, where $k \in [-0.5, 0.5]$

then $\phi$ for the strongest bin 10 looks like $k·\pi$.

Some sample measurements:

$k=-0.5 \Rightarrow \phi=-0.499·\pi$

$k=-0.25 \Rightarrow \phi=-0.252·\pi$

$k=-0.15 \Rightarrow \phi=-0.151·\pi$

$k=-0.05 \Rightarrow \phi=-0.054·\pi$

The deviation gets bigger as you move towards the bin, though it's pretty good near the middle. Is this some kind of a parabolic correlation?

So what's the math behind this? Why do these values seem correlated?

And if they are correlated, could this be used for finding the peak frequency with some precision?

For a set of exact equations that describe how a pure real tone works in a DFT check out equations (23)-(25) in my blog article:

DFT Bin Value Formulas for Pure Real Tones

The phenomenon you have found can be seen clearer if you consider the complex definition of a real valued sinusiodal.

$$ \cos( \theta ) = \frac{e^{i\theta}+e^{-i\theta}}{2} $$

So, when you are near a peak bin, one of the numerator values is dominant, and you have an approximation of a complex tone.  Using equation (24) in this article shows how a the distance to the bin $\delta$ rotates the bin value.

DFT Bin Value Formulas for Pure Complex Tones

Hope this helps.

To answer your last question: Yes, I've done it different ways with more coming.

Check out:

Exact Frequency Formula for a Pure Real Tone in a DFT
Two Bin Exact Frequency Formulas for a Pure Real Tone in a DFT
Improved Three Bin Exact Frequency Formula for a Pure Real Tone in a DFT
A Two Bin Exact Frequency Formula for a Pure Complex Tone in a DFT
Three Bin Exact Frequency Formulas for a Pure Complex Tone in a DFT
[Jacobsen's Frequency Formula as an Approximation]
[Candan's Tweak of Jacobsen's Frequency Formula]
[Tweaking Macleod's Frequency Formula]

Your $k$ is the fraction of a bin distance from a bin on the frequency scale of cycles per frame.

Thus, from (20):

$$ k = -\frac{N}{2\pi} \delta $$

The negative sign is unfortunate, I should have done it differently in the article.

$$ \delta = -\frac{2\pi}{N} k $$

Your $\phi$ is the angle of the nearest complex bin value, so I will use $\theta$ in the formula for the phase angle from the signal definition in (24):

$$ \phi = \arg(X[]) = -\delta (N-1)/2 + \theta $$

$$ \phi = \frac{2\pi}{N} k (N-1)/2 + \theta $$

$$ \phi = \frac{N-1}{N} k \pi + \theta $$

So, for N being large, and $\theta$ being zero, $\phi$ looks like $k \pi$.

Remember, this is based on using a complex tone as an approximation for the real value tone, so it won't be exact.  I am not sure why your ratios are slightly larger than one when they should be slightly smaller.


MATLAB appears to be applying my filter to the spectrum twice, once on each side So it looks like MATLAB is applying the filter to the "negative" and "positive" frequencies separately and I end up with two band pass regions with a null around zero. I want the whole spectrum filtered once.

[b,a] = butter(2,[.18 .3739],'bandpass');
hm = filter(b,a,Y);
figure;
plot(20*log10(abs((fftshift(fft(hm))))))

Any ideas?

Your filter has a real-valued impulse response, which means it has a two-sided spectrum.  Your filter is only applied once, but that application changes both positive and negative frequencies: -0.3739 to -0.18 and from 0.18 to 0.3729 (modulo the transition region of the filter).  If the filter did not do this, you would have a low pass filter and not a bandpass filter.


Find the error in phase angle and dB gain Function Transfer This is a question paper that I've tried to solved but seems I am missing something.

My attempt:

$$ |G|_{w=0} = K$$

So,

$\require{cancel}|G|_{w=0.5a} = \dfrac{K}{\left|1+\dfrac{j0.5\cdot \cancel{a}}{\cancel{a}}\right|} = \dfrac{K}{\underbrace{|1+j0.5|}_{Gain}}$

$Gain = 20 \log(\sqrt{1,25}) = 0.96910013008$ [OK]

But about the phase error?

$\angle(1+j0.5) = \tan^{-1}(0.5/1) = 26.56º$

EDIT: Sorry picture this question:

What I did wrong?

The "Error" would be the actual phase and amplitude relative to the asymptotic plot shown. You computed the gain and phase, so now compare this to the asymptotic lines given by the plot to determine the error.


What is the procedure to obtain a 3D face dataset? Can anyone suggest a good 3D face dataset? Project Description:

Input:

Frontal Face Image
Expression
Angle
Size

Details:

I have to convert the input frontal face image into 3D face and simulate the given expression (smile, cry, neutral, ..); rotate by the given angle and re-size to the given size. Again convert back into 2D face after the changes.

After googling, I found that Morphable models is a good algorithm to start with; I want a 3D face database to start & complete the implementation.

I found some 3D databases list on www.face-rec.org. Can anyone help me with buying (pricing, procedure ....) the database, and how to choose the database?

Thank you in advance..

I don't personally know of a good one to vouch for, but on face-rec.org/databases/
at least, isn't the "GavabDB: face database" free?

see:
http://gavab.escet.urjc.es/recursos_en.html

has 3 (apparently free) downloads:
Brief description: GavabDB is a 3D face database. It contains 549 three-dimensional images of facial surfaces. These meshes correspond to 61 different individuals (45 male and 16 female) having 9 images for each person. The total of the individuals are Caucasian and their age is between 18 and 40 years old. Each image is given by a mesh of connected 3D points of the facial surface without texture. The database provides systematic variations with respect to the pose and the facial expression. In particular, the 9 images corresponding to each individual are: 2 frontal views with neutral expression, 2 x-rotated views (±30º, looking up and looking down respectively) with neutral expression, 2 y-rotated views (±90º, left and right profiles respectively) with neutral expression and 3 frontal gesture images (laugh, smile and a random gesture chosen by the user, respectively).

Note:
Those publications that use this signature date must reference the following work: A.B. Moreno y A.Sanchez. GavabDB: A 3D Face Database. Proc. 2nd COST Workshop on Biometrics on the Internet: Fundamentals, Advances and Applications, C. Garcia et al (eds): Proc. 2nd COST Workshop on Biometrics on the Internet: Fundamentals, Advances and Applications, Ed. Univ. Vigo, pp. 77-82, 2004

Readme.rtf (9 kB; .rtf)

CD1_faces1-20.zip (38 MB; .zip)

CD2_faces21-40.zip (39 MB; .zip)

CD3_faces41-61.zip (38 MB; .zip)


Phase diagram of a rectangular pulse with Fourier Series - help understanding I understand perfectly fine how to plot the magnitude of a Fourier series, but I'm having serious trouble understanding how to plot the phase spectrum. Below is a picture of a rectangular pulse. The magnitude diagram makes sense to me - but how did they determine the phase diagram?

I know that the coefficient is:
$a = \frac{Asin(n*pi/2)}{n(pi/2)}$ or $Asinc(n*pi/2)$, but I'm not sure how to determine the phase spectrum from this.

I thought of using Euler's formula to convert the sine function, but then I'm a bit lost on how to use the exponential values to find the phase (if that's even correct?). So I got it to look like this:
$$a = \frac{A(e^{jn*pi/2}-e^{-jn*pi/2})}{(2jn*pi/2)}$$

But then.. I'm not really sure? Because in the diagram the angle is pi, but I can't see how to get that.

The phase of a complex number is the angle that the number makes with the positive real axis on a complex(Imaginary vs real) plot.

Since the coefficients are given by $$a=A \frac{\sin(n\pi/2)}{n\pi/2}$$ it can be seen that it is a real number for any value of n. So $a$ will have a phase of either 0 (for positive number) or $\pi$ (for negative number). Hence whenever the value of $a$ becomes negative the phase is $\pi$.


Fast Fourier Transform showing the same results but with different files? so I'm pretty sure that I'm plotting the wrong frequency domain against my FFT values because I keep getting basically the same graph shape every single time. I think that I just don't comprehend the context of some of the values, as I'm quite new to signal processing. Attached below is the frequency domain I'm attempting to use, which is what I'm seeing on most websites. Is T the actual period of your data? It doesn't really matter whether or not I have user input or hardcode a value in, the same graph shape remains. Same goes with completely different data sets, always the same shape remains.

import matplotlib.pyplot as plt
from scipy.fft import fft, ifft

transform = fft(data)
per = input("What's your period: ")
T = int(per)
f = np.linspace(0, 1 / T, N)
plt.plot(f, np.abs(transform)*(1/N))

I'm not familiar with this type of code, but I can take a guess at what you want to do. Your variable names are inappropriate. I suggest you change  your code to something like the following (comments are inside brackets):

spectrum = fft(data)   [Assumes you've transformed a data sequence whose length is 'N' samples.]

fs = input("What's the 'data' sampling rate (in Hz): ")

delta_freq = fs/N   ['delta_freq' is the frequency spacing (measured in Hz)
between the spectral samples contained in the 'spectrum' sequence]

freq_axis = np.linspace(0, delta_freq, (N-1)*delta_freq)

plt.plot(freq_axis[:N//2], np.abs(spectrum)[:N//2]*(1/N))  [Plots the normalized (i.e., the 1/N scaling) spectral magnitude over the positive-frequency range]


Combining two filters (audio processing) I'm trying to combine two filters into a single filter command for a homework problem and was looking for advice. I don't necessarily need an explicit answer as long as someone can point me in the right direction.

In the following code, assume that 'y' is a sound file that has been pre-loaded (10kHz). Also, both the hlpf and hnotch filters work together separately to produce an audio signal with noise and a beeping tone removed. However, when I try to combine the two into a single filter, I do not achieve this result. Here is my code:

Fs = 10000;
f0 = 1000;                %#notch frequency
fn = Fs/2;              %#Nyquist frequency
freqRatio = f0/fn;      %#ratio of notch freq. to Nyquist freq.
notchWidth = 0.1;       %#width of the notch

%Compute zeros
zeros = [exp( sqrt(-1)*pi*freqRatio ), exp( -sqrt(-1)*pi*freqRatio )];

%Compute poles
poles = (1-notchWidth) * zeros;
b = poly( zeros ); %# Get moving average filter coefficients
a = poly( poles ); %# Get autoregressive filter coefficients

hnotch = [poly(poles)];
alpha = 0.42;
hlpf = fir1(100, alpha);
hcombo = conv(hlpf,hnotch);
s2 = filter(hcombo,1,y);
soundsc(s2,Fs)

Any help would be greatly appreciated.

In order to combine two filters, you need to convolve their impulse responses, and you haven't done that in this case.

In fact, according to this code snippet, the notch filter is incorrectly built.

I'm sure once you review it you'll figure out your errors. Good luck!


How do I use a Savitzky Golay filter to find local maxima (in between samples) in a discretely sampled 1D signal? I have a seismic signal y(i):

Here I have found one maximum: i=152.54, y=222.29 manually and plotted it in red.

I want to find all maxima automatically.

I read that the Savitzky Golay Filter (SGF) can be used to find smoothed estimates of both a signal and its derivatives, and that one of the benefits of the SGF is that it preserves minima and maxima much better than other filters. This sounds great for my use.

I found a Matlab script that generates SGF coefficients.
And used this to find that the 4th order SGF coefficients for the derivative.
I coded a small Matlab script that

finds the derivative of the signal by convolving the signal with the
4th order SGF coefficients for the derivative
finds pair of samples (i,i+1) where the derivative changes sign
finds zero crossing of derivative by linear interpolation between i and i+1

Script:

function [maxX,maxY] = findLocalMax(y)
% Kernel for 4th order Savitzky-Golay filter for finding derivative:
d4 = [0.0724 -0.1195 -0.1625 -0.1061 0 0.1061 0.1625 0.1195 -0.0724];

dy = conv(y,d4,'same'); % derivative

[m n] = size(dy);
maxX = [];
maxY = [];
for i = 1 : n - 1
  if dy(i) < 0 && dy(i+1) > 0 % max somewhere between i and i+1
    a = dy(i)/(dy(i) - dy(i+1)); % linear interpolation
    mx = i + a;
    maxX = [maxX mx];
    my = y(i)*(1-a) + y(i+1)*a; % linear interpolation
    maxY = [maxY my];
  end
end

Im my script I had to test for if the derivative changes from negative to positive to get the function to give the desired result, however this confuses me.
Should not the derivative for a maximum go from positive to negative?
Is there some better way to distinguish between maxima and minima?

Below is the result of using this function to find the maxima on my signal:

The results look good, but I notice that some maxima are not found:
i= 143.13, 190.88, 256.97.

Is this because they are to close to other maxima?

How can I control the nearest two maxima may be?

Thanks in advance for any answers!

While I'm not familiar with this specific type of filter, based on the plot you've shown, I would guess that the maxima that aren't found by your process are just butting up against the time resolution inherent in the process. Any kind of "smoothing" implies that there is some time-local smearing of the signal of interest, such that if there are two nearby peaks, it's possible that they will be coalesced into one. It's possible that a lower-order filter might present less of this behavior, likely at the expense of the amount of smoothing that you get.


Control optimization problem I am running into a problem where I have a control system $S[t]$ that takes a control $C[t]$, so that
$$S[t+1] = H(C[t,t-1,...], S[t,t-1,...])$$ the response of the system is the history of controls and responses.

I am looking for a way to control $C[t]$ to optimize the value of $S[t]$ at any point in time.

What is the literature on this problem called? Where do I start looking? Thanks.

This looks like it might be able to be written as a dynamic programming problem. Bellman's equation looks very similar to the way you have expressed your control problem.


the importance of discrete sine transform what are the applications of discrete sine transform . can discrete sine transform be used in place of discrete cosine transform ? please state the reason.

can it be used in compression.?

regards,
phani tej

DCT performs close to KLT so that it is the best choice to compress highly corelated signals. For low correlated signals DST is used , check video coding standards HEVC. Therefore, one is the complement of the other, never a replacement. DCT has the DC component and DST has no DC component


What Is an Oriented Gaussian Second Derivative Filter In the paper: Detecting and Localizing Edges Composed of Steps, Peaks and Roofs, the authors refer to an image filter as an oriented second-derivative Gaussian filter. I'm trying to figure out what this means.
From my understanding a Gaussian filter for a given standard deviation $\sigma$ of size $n \times n$ for some odd positive integer $n$ is given by the formula $g(x,y)  = \frac {1}{2 \pi \sigma^2} e^{\frac{-(x^2+y^2)}{2 \sigma ^2}}$ applied to image co-ordinates of a rectangular region of a 2-D image (ignoring padding for now).
Then there are Guassian first-derivative filters consisting of $g_x(x,y)$ and $g_y(x,y)$, the partial derivatives, which for a given angle $\theta$ can be combined into an overall oriented filter, say:
$t(x,y, \theta) = cos (\theta) g_x(x,y) + sin (\theta) g_y(x,y)$
as per the answer to my question here.
Now it seems to me there are some choices for what could be considered by the term oriented second-derivative Gaussian filter (which after some Google searching I could not find a definition of):
a) An orietned Laplacian of Gaussian (since it involves second derivatives), a.k.a $L(x,y, \theta) = cos (\theta) g_{xx} (x,y) + sin (\theta) g_{yy}(x,y)$ (where $g_{xx},g_{yy}$ are the partial derivatives twice with respect to $x,y$ respectively). This seems the most likely choice.
b) Some kind of mixed partials derivative filter like $t(x,y, \theta) = cos (\theta) g_{xy} + sin(\theta) g_{yx}$, although the mixed partials should be equal by Clairaut's theorem.
Any insights appreciated.

Unless mentioned otherwise withing the context the classic interpretation of Second Derivative Gaussian Filter is indeed (a) in your question:
$$ L \left( x, y, \theta \right) = \cos \left( \theta \right) {g}_{xx} \left( x, y \right) + \sin \left( \theta \right) {g}_{yy} \left( x, y \right) $$


Solving LASSO (Basis Pursuit Denoising Form) with LARS I'm now working on using LARS (Least Angle Regression) algorithm to solve a LASSO problem in Basis Pursuit Denoising form like:

\begin{align*}
\quad && \arg \min_{\beta}{\left\| y - X\beta \right\|}_{2}^{2} + \lambda {\left\| \beta \right\|}_{1} && (1)
\end{align*}

I have done many search jobs on the web, what confused me is LARS is more like an algorithm to solve problems in this form:

\begin{align*}
\quad & \arg \min_{\beta} {\left\| y - X\beta \right\|}_{2}^{2}
&& \text{subject to.} && {\left\| \beta\right\|}_{1} \leq t && (2)
\end{align*}

I know these two forms is equal to each other in theory and LARS can solve problem(1) for all 0<λ≤∞.
But I want know how to use LARS to solve problem(1) in practice, which means how can I write a function like:

function lars(X, y, lambda){
 ...
 return beta;
}

There are 2 forms of the Basis Pursuit problem:

$$\begin{align*}
\text{The $ \lambda $ Form:} & \quad && \arg \min_{x} \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2} + \lambda {\left\| x \right\|}_{1} \\
\text{The $ \epsilon$ Form:} && \quad & \arg \min_{x} {\left\| x \right\|}_{1} \\
& \text{subject to} &&  \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2} \leq \epsilon
\end{align*}$$

Real world model is:

$$ A x = y $$

Where $ x $ is a sparse vector.
Yet, in reality we don't have measurements of $ y $ but $ b = y + v $ where $ v $ is a vector with the properties of our measurement method.

Hence we allow the model not to have strict equality which implies:

$$ {\left\| A x - b \right\|}_{2}^{2} \leq \epsilon = {\left\| v \right\|}_{2}^{2} $$

Now, the different models are equivalent as for any $ \epsilon $ there is a $ \lambda $ (Which depends on $ A $ and $ b $ unfortunately) which the models ( (1) and (2) ) are equivalent.

For instance I created simple simulation on for that simulation:

The full code is available on my StackExchange Cross Validated Q291962 GitHub Repository.

Since you linked to the paper you're reading - Sparse Geometric Representation Through Local Shape Probing the problem solved there is the $ \lambda $ form.
Since it is a strict Convex Problem there is a single solution. LARS isn't considered a very efficient or fast method to solve the problem, hence I recommend using a different solution.

I have a project with many solvers of the $ {L}_{1} $ Regularized Least Squares. I suggest you just take the Proximal Gradient Method - SolveLsL1Prox.m or the Accelerated Proximal Gradient Descent Method - SolveLsL1ProxAccel.m.

So you can just pick any of those MATLAB codes and solve the problem.


3D reconstruction of an image based on depth of intensity levels? Is there anyway of segmenting objects based on how far they are from the viewer?

Can color values assess this sort of a thing for us?
How will intensity levels help in determining how far the object is from the viewer?

Another image:

This is an extremely difficult task, one which is a very active line of research. I've managed to find a semi-recent paper on the subject. I won't go into the details, but here's a few things that this paper found that can be used.

Textures: If you see a given texture, they can help you figure out how far something is away. The wood texture here would be quite ideal for this.
Fixed object sizes: Identify common objects, see how large they are, and use them to help with depth analysis. The phone book, soda bottle, etc should work well with this.
Focus: If you know how the picture was take, focus can help you identify ranges.
Shadows- If you know where the light was coming from, shadows can help you identify how tall objects are, and thus how far they are away.


How does sensor fusion work? I want to understand the math/physics/algorithms Sensor fusion algorithms can give a more precise 3D orientation (and possibly postion?) of a device by combining readings from an accelerometer, gyroscope, and magnetometer.

Can anyone explain, or provide links, explaining the details behind this? I want to understand the math and physics, such that if I have raw 9-DOF sensor data, I might be able to implement sensor fusion. Assume I have sufficient background of linear algebra, calculus, etc.

Cheers!

Sensors provide analog or digital insights into a reality that is difficult to grasp. Just like your senses do.

Sensors are designed to seize a special portion of a measurable phenomenon: electrical, chemical, physical...  Just like your senses do. The usual five ones are sight, hearing, touch, smell and taste, but there are more to human skills.

Understanding a physical phenomenon via signal processing is like trying to tell the main ingredients from the  recipe of the meal you ordered in a restaurant last evening. You do not know what happened in the kitchen, yet you have chosen what to order, and have seen, tasted, smelled, even touched the plate.

Each of your senses and understanding provided you with hints or knowledge, but you can never be sure, as some cooks master sense deception. What you see as an artifial green marmelade made of strawberries can taste like apple jelly.

Based on you cooking expertise, the combination of all your senses can bring you closer to the the actual recipe, provided that:

you can use your senses often enough: if you take only one bite every 10 days, it is unlikely you get the menu. This is data sampling;
each sense is sensitive enough. This is sensor sensitivity;
your sense list is close to complete, and the coverage is sufficient for your purpose. If you cannot taste sugar anymore, you won't be able to cook some meals for people who still can taste it. This is measurement span;
you are able infer, or model some of the processes you are unaware of. This is modeling.

Sensor fusion is the art and science of combining sensory data, knowledge and models from disparate sources so that the resulting information has more validity or less uncertainty than individual  sources.

Examples: The act of averaging $N$ sensor data for the same deterministic signal, with different stochastic realizations of a noise term, is the most basic sensor fusion operation. With the model of independent Gaussian noises of deviation $\sigma$, the average theoretically yields a $\frac{\sigma}{\sqrt{N}}$ deviation(less uncertain). The combination of three channels (Red, Blue, Green) offers a color image (more valid than a monochannel gray-scale).

The first example combines information from a single sensor model. The second one from sensors operating in the same domain (electromagnetic waves), yet in different portions of the spectrum. In general, fusion operates on different sensors with different rates, ranges, domains and mostly units.

The problem is highly dependent on the phenomenon you are looking at, the available sensors, and the information you are looking for.

The physics tell you the potential information you can get from your sensors. The mathematics can model how they are related or complementary, or what is irrelevant information (noise). The algorithms will combine the previous knowledge as optimally as possible, in terms of precision, accuracy or speed.

The topic is related to the realms of Sensor fusion, Data fusion or Information integration, with a short overview in Principles and Techniques for Sensor Data Fusion. Many more books are available, as in Best book for learning sensor fusion, specifically regarding IMU and GPS integration.

Regarding your actual problem, a first step would consist in understanding What are the differences between a gyroscope, accelerometer and magnetometer? which could help you push techniques a little further. And achieve the goal of fusion: using sensor differences, sum their data, in the most clever way.


numpy.correlate and autocorrelation; audio signal I want to calculate delay between an input and an output audio signal of my audio processing system. The input and output signals are available as signed 16 bit integers. To try out, I tried the autocorrelation of the input signal with the following numpy commands:

import numpy as np
import wave
wfp = wave.open('test.wav', 'rb')
samples = wfp.readframes(wfp.getnframes())
signal = np.frombuffer(samples, np.int16)
corr = np.correlate(signal, signal, "full")

I assumed that the peak of the autocorrelation is always at lag 0, which is at the index corr.size / 2 in the corr array. However I get different values for the index, when I calculate np.argmax(corr).

When I normalize the signal first to get values between -1.0 and +1.0, the peak is always at corr.size / 2 as expected, at least in my tests with different signals.

For normalization, I used the following steps:

signal = signal / float(0xFFFF)

Can someone please explain to me, why I have to normalize the signal.

Having had the same problem and without success to find a tool to sync the start of video/audio recordings automatically,
I decided to make syncstart (github).
The basic code is this:
import numpy as np
from scipy import fft
from scipy.io import wavfile
r1,s1 = wavfile.read(in1)
r2,s2 = wavfile.read(in2)
assert r1==r2, "syncstart normalizes using ffmpeg"
fs = r1
ls1 = len(s1)
ls2 = len(s2)
padsize = ls1+ls2+1
padsize = 2**(int(np.log(padsize)/np.log(2))+1)
s1pad = np.zeros(padsize)
s1pad[:ls1] = s1
s2pad = np.zeros(padsize)
s2pad[:ls2] = s2
corr = fft.ifft(fft.fft(s1pad)*np.conj(fft.fft(s2pad)))
ca = np.absolute(corr)
xmax = np.argmax(ca)
if xmax > padsize // 2:
    file,offset = in2,(padsize-xmax)/fs
else:
    file,offset = in1,xmax/fs


Subcarrier assignment in MIMO-OFDMA system Is it possible for multiple users to simultaneously transmit on the same subcarrier in a MIMO-OFDMA system?

For a given OFDM symbol, no subcarrier can be assigned to two users simultaneously. Each subcarrier is uniquely assigned to a unique user. Else you have to implement an Interference cancellation at receiver (which would require having channel knowledge of other users as well, meaning global CSI) or at the transmitter and extremely complex interference alignment/precoding has to be done.

In practical systems like LTE or 5G, thus one subcarrier in a given OFDM symbol is allocated to only one user for the reasons listed above.


How to choose basis functions that contribute most efficiently per term to an approximation of an image f(x,y)? GOAL:

I would like to approximate some positive, scalar function, $f(x,y) > 0$, on a 2D field of finite size i.e. $x=[a,b],y=[c,d]$

OBSTACLE:

I am familiar with the set of basis functions used in the Fourier series but using the Fourier series requires too many terms. Infinitely many:)

Instead, I desire to find a collection of functions that can be used to approximate $f(x,y)$ with a "finite" number of terms.

QUESTION:

How can I approximate $f(x,y)$ with a collection of functions $G(x)$ and $H(y)$ (see definition of $G$ and $H$ below), such that the number of elements in $G$ or $H$ is finite?

$$G(x) = [g_1(x)\ldots g_n(x)]\\
H(y) = [h_1(y)\ldots h_n(y)]$$

The elements in $G$ and $H$ do not have to be of the same form ($g_1(x)=\sin(nx),g_2(x)=\sin(mx)$)

Thanks in advance:)

EDIT:

I should clarify that $G(x)$ and $H(y)$ need to be combined in an inner product fashion i.e.

$$f(x,y) = < G(x),H(y) > = g_1(x)h_1(y) + g_2(x)h_2(y) + ... g_n(x)h_n(y)$$

I do not  understand why you need an infinite number of  terms for a Fourier approxiimation. Yet, if you are looking for tools on non-discretized 2D data,  I would suggest you to start with the concept of Matching pursuit of images. You take a dictionary of 2D atoms that is complete (could be complex exponentials, gaussian bumps). Then, you find the  atom that best matches you  data, subtract it with the necessary scale factor, and repeat.

You could also resort to anisotropic  triangulations, but  this would not suit your "need to be combined in an inner product fashion".

The  quality of approximation could be evaluated over some classes for regular functions ($C^p$ smooth inside $C^q$ curves), and are related to properties of  your  dictionary. You might find some references in section 4.1.1. Matching pursuits from A panorama on multiscale geometric representations, Signal Processing, 2011, L. Jacques et al.

I do suspect that the inner-product  constraint will limit the efficiency: optimal approximations in 1D do not extend generically  in 2D.


What is "Equivalent Time Sampling" and what is it good/used for? The literature on this method seems scarce, but I know that it has been used on radar systems to 'get away with' not having to sample nearly as fast as Nyquist would otherwise dictate, (at the cost of additional time). Or so I think.

So, how does "Equivalent Time Sampling" work, and where can it be used? How does it square with the Nyquist Criterion?

Note that the Nyquist criterion applies to bandlimited signals; but the so-limited band does not need to be baseband, or even contiguous, just not overlapping (aliased) with any other sub-bands after folding.

Periodic sampling fan-folds an infinite spectrum in to "fans" with a width of half the sampling rate.  Any frequency bands, after this fan-folding, that do not overlap with any other frequency bands in the spectrum when so folded can be reconstructed, according to one form of Nyquist theorem.

In practice, any random phase jitter in undersampling will cause the bands too far above Fs/2 to become more phase noise than signal.  If the phase noise is mostly white and uncorrelated with the signal, then sampling for a much longer time frame may allow averaging out this phase noise sufficient to reconstruct a repetitive narrow-band signal, with some probability of accuracy for a given extended sampling time window.

So, some oscilloscopes have a mode where they undersample a bandpass-filtered narrowband input, then store and composite a much longer sampling time by folding the samples down to baseband.

In addition, signal reconstruction of a bandlimited signal may be possible with non-periodic sampling of known sample times, as long as the local density of samples is high enough.  This is numerically much less stable, as estimating a reconstruction of the spectrum or waveform from a minimum number randomly spaced samples becomes more of a sensitive linear optimization problem, rather than a closed form transform or convolution.  However some oscilloscopes might use a pseudo-random non-periodic undersampling, but averaged over a long time period, and let the user determine when the density of samples "looks good enough" when reconstructing an image of an undersampled repetitive waveform.

Randomized sampling can avoid certain aliasing problems, such as perfectly periodic sampling not allowing reconstruction of an arbitrary phase sinusoid with a frequency exactly at Fs/2 (it's smooshed at the creases in the fan-fold).


Correlation filter output range normalization I'm developing correlation filters based image recognition.
I implemented MACE correlation filter in matlab:

training code:

D = diag(mean(abs(X),2));
% inv(A) * B = A \ B
XDX = ctranspose(X) * (D \ X);
h = (D \ X) * (XDX \ u);
H = reshape(h, size(I));

where X is the d.^2xN matrix of vectorized FFT2 transformed training images (images FFT2 transform spectra are vectorized by concatenating columns), u is the vector of image membership (1 for images to recognize and 0 for the rest of the samples), I is the dxd sample training image to reshape correlation filter h to H matrix.

testing code:

R = ifftshift(ifft2(I.*conj(H)));

where I is the dxd test image zero mean and unit variance normalized and H is the trained correlation filter.

I trained one correlation filter on 61 training images from which there is just 1 image to recognize with 1 in vector u and remaining images with 0's in vector u.

However despite reported results in the papers on correlation filters where correlation plane output is in range 0 ... 1 my results are up to 10e-4.

Is there any normalization typically applied to scale correlation plane output?

Besides I found that without conj(H) in the testing code there is no correlation peak for images to be recognized. I was unable to find in the papers that missing operator in the formula which I found at mathworks forum

I was just working on the same issue and the magnitude was the same as yours (10e-4), but the peak in the middle, in terms of the shape, looked correct.

I was able to solve this problem by scaling back R, by d, the number of pixels in the image.  In my case I am using 64 x 64 pixel images (d = 64*64 = 4096).

Also, make sure you take the real of R, to get rid of the complex part, even though it should be very small, then multiply R by d and your scaling should be correct.  To test this, use an image you trained your filter on and your max(R(:)*d) = 1.0.

% COMPUTE THE CORRELATIN PLANE
R = ifftshift(ifft2(test_img.*conj(H_square)));
R = real(R);
% PLOT THE CORRELATION PLANE
figure;
mesh(R*d)
axis([0 64 0 64 -0.5 1]);


Why not use synthetic datasets for training machine vision (deep learning) algorithms? I am studying some examples of how to perform 3D instance segmentation in indoor scenes, and I have noticed many of the available datasets are from real environments. I was wondering with the availability of animation and computer graphics why people don't use these methods to create datasets ?
For one thing they are much, much simpler to annotate (for example for instance segmentation) since we already know the objects in the scenes. And we will have a lot more variety in the dataset since we have a much bigger freedom in in choosing and manipulating the objects in the scene (type of objects, color, lighting).
I do understand that for some more advanced tasks like analysing how humans interact with the environment the animations would be less realistic, but for tasks like instance segmentation it seems to me that the datasets would be accurate enough.
I do understand that creating these scenes (specially if we want to create an animation) could also be time consuming but wouldn't that be offset by the ease of annotation ? (specially for 3D scenes)
So what am I missing here ? why don't we have more of these types of datasets ? what is the drawback of using such a dataset ?

Well, when you're training something against a dataset, you're training it to work on that kind of data.
Since we have no idea of what the neural network "uses" from the input images to fulfills its purpose, these might be things where our computer graphics aren't quite like real-world images.
We're "trained" with real-world imaging, perceived through our eyes with our brains. The neural net has none of that. It might simply do the segmenting, and that's actually pretty likely, on the specific kind of anti-aliasing done by the 3D renderer. Great, now you have a neural net that can segment generated images very well, but can't deal with real-world photography. Not quite sure what you've gained there!
But: sometimes, you actually do that, if you can describe your real-world data well enough – but that's usually things like RF signals, abstracts datasets etc., and often you'd intentionally add noise, so that the noise equivocation is more prominent than the differences between synthesized data and the real world.
You'd always have to verify, either way.


How i can smooth a sinusoidal signal using it's local maximum I have a signal that it is like a Sinusoidal signal with many local maximum, I want to smooth this signal with almost connecting it's local maximum to each other.

You could use the Hilbert transform to calculate the signal envelope. It is discussed in this question.


General tips for PID tuning of super low-friction actuator? I realize this question is not directly related to signal processing, however, it's relevant to system analysis which is relevant to most signal processing engineers. There's also no good alternative stack exchange for me to post this question to (please let me know if there is though).

I'm currently struggling to tune an air-bearing actuator. Because this actuator has such low friction, most of the energy put into it will not dampen itself out naturally. Instead, I'm relying on the derivative gain of my PID controller to do that.

However, I'm struggling to achieve a stiff system (high proportional gain) without it going unstable. I'm using a high powered amp with this stage which is plenty capable of creating stiff control, however, when ever I increase the proportional gain the stage goes unstable. My hope was that I could increase the derivative gain in response to this, yet even doing this it doesn't seem to make the system stable.

What I'm running into is a system that is "squishy" controlled that then goes unstable when I try to increase the proportional gain.

What exactly am I doing wrong?

1 - Do you have some general idea of what your transfer function looks like? Or at least a frequency response? Could you perform some kind of step response? You could then try to fit an order-2 model on your step response and then try to tune your controller.

2 - You mention that your system is underdamped. I assume that you have underdamped stable poles (in the left-hand plane) in your process transfer function. Perhaps you could place zeroes near the poles to damp the system.

3 - I never use a PID with a "real derivative" i.e. derivative = x[n] - x[n-1], the derivative is always filtered in order to stabilize the system. Otherwise modeling errors and high-frequency noise will be amplified by the derivative which is something you don't want. Or instead of a PID you could use a PI with a lead controller cascaded.

4 - I recommend that you consult this book http://www.cds.caltech.edu/~murray/books/AM05/pdf/am08-complete_22Feb09.pdf