Compressive Sensing with Square Measurement Matrix If my measurement matrix have same number of row and column and the unknown vector is sparse can I still use Compressive Sensing to get better reconstruction with fewer measurement?

It's important to clearly distinguish the terms Compressive Sensing (CS) and Sparse Signal Recovery (SSR). CS is about taking fewer measurements than what classical criteria such as Nyquist would dictate us. In the discrete setting this means taking $M<N$ measurements of an $N$-sparse vector.

SSR is about how to recover our $N$-sparse vector from your measurements.

That said, this should answer your question: If you apply CS, then you need SSR to recover your signal, however, the converse is not true - you can apply SSR even if you didn't do any subsampling. For SSR, your measurement matrix can be flat or square, it could even be tall.

Mathematically speaking, for SSR you solve a problem of the sort $$\min_x \|y - A x \| + \lambda h(x),$$ where $h(x)$ is a regularizer that encourages sparse solutions, i.e., an $\ell_0$-"norm" or its relaxation, the $\ell_1$-norm. In CS, such problems are solved since $A$ is flat and hence we must regularize to find our sparse solution out of the infinitely many possible solutions. The same problem can be used in a non-CS setup where $A$ is square or even tall, e..g, to find an approximate solution that shows a desired level of sparsity.


How to use Principal component Analysis(PCA) to extract feature from images How can use principal component Analysis to extract feature from image that is been acquired by the camera?

PCA is mostly used to reduce the dimensionality of the features.

You should first use a feature detector before you can use PCA to reduce the dimensions of these features. Common feature detectors are SIFT, SURF, ORB etc. Some of them are even available in OpenCV.


Use Recursive Least Square as an "Actual" Filter? I am currently studying recursive least square (RLS) and as far as I understand, the setup is (for example) that we are given a process with uncertain parameters/noise and we want to adapt parameters of a model of that process online.

To me, that sounds like a Luenberger observer, where the parameters of the underlying model are adapted online, depending on the model quality.

If this is true, a side question would be why this technique isn't usually called an observer, but a filter instead? However, my main question is:

Question: Can RLS also be used instead of a standard low-pass filter?

So, given only a signal from which I want to remove high-frequency components (like noise), can a RLS filter be used to do that as well? And if so, would it have any advantage over normal low-pass filters?

The problem that recursive least squares (RLS) can solve can be formulated as recursively solving for $\hat{\theta}$, such that it is the least squares solution to

$$
\hat{\theta}_n = \arg\min_x \sum_{k=0}^n w[k]\,\|z[k] - \phi[k]^\top x\|_2^2,
$$

where $w[k]$ are weights, $z[k]$ and $\phi[k]$ are known and $z[k]$ is assumed to be generated by using $\phi[k]^\top \theta$ with $\theta$ constant and unknown to the user. The recursive part of RLS refers to that the solution $\hat{\theta}_n$, plus additional variables that encode the problem, are used to calculate the solution $\hat{\theta}_{n+1}$. Now, if $\phi[k]$ satisfies a persistence of excitation condition, then over time $\hat{\theta}$ can be shown to converge to the value of $\theta$. A more detailed formulation and derivation can be seen in Mathematics StackExchnage - What's the Idea Behind Recursive Least Squares (RLS)? How Could It Be Derived?

So this formulation initially suggests that it would only be suited for estimating constant variables. However, by using exponentially decaying weights (so weigh old data exponentially less and less) RLS can still behave similar to a low-pass filter. So that the RLS can still reasonably track a slowly varying $\theta$.

It can also be shown that RLS (with constant weights) is equivalent to a Kalman filter for the system with $F_k = I$, $B_k = 0$, $H_k = \phi[k]^\top$, $Q_k=0$ and $R_k=1$. So if you want to filter out noise from a time varying signal and you have a model of how you think that signal should change over time, you might be better of using a Kalman filter using that knowledge. Namely, RLS with exponentially decaying weights will always lag behind a changing signal, similar to a low-pass filter, while a Kalman filter corrects for this using the prediction step.

Finally, regarding RLS being called a filter and not an observer might be because the label observer is I believe usually reserved for estimating properties of dynamical processes, which the constant $\theta$ which RLS tries to estimate is lacking. However, I am often also confused about when to use the term filter, estimator or observer, so difference in meaning might be small.


Image processing free open source libraries This is the first time that I post a question  here so excuse me for asking such a beginner question.
what free open source 3rd party libraries which have good documentation(because I'm new to using 3rd party libraries) would you suggest me for fourier transform,inverse fourier transform , meshgrids , dynamic arrays and etc to have a better performance and higher speed in visualC++ when implementing image processing filters?
Note that I'm not permitted to use openCV. In fact I'm going to compare the results of my code with results of openCV

For Fourier transforms go get fftw.  I don't know what meshgrids are, so I can't answer that.  For dynamic arrays: that's what std::vector is, it's built into the language.

For a wider variety of image transformations you might look at GEGL and the library underlying Imagemagick.

There's also this question about filtering software: Use Butterworth and Chebychev filters.


Sharpening an image using emgu I have c# desktop application and I am using emgu framework as a wrapper for OpenCV.

I am trying to sharpen an image because the input images are not great.

I tired a technique I found via searches that suggest from the original image:

I take a copy and blur it.
I then take another copy of the original and invert it.
I then merge the 2 new images together.
I then finally invert that merged image.

I have done all that but the image pretty much looks the same.

Is this the right approach?

I can post my code if need be - I used Aforge to quickly test this approach though.

UPDATE:
To answer my own question I used this:

            Image<Bgr, byte> image = new Image<Bgr, byte>(@"D:\20140320022038047.jpg");
            Image<Bgr, byte> image2 = new Image<Bgr, byte>(@"D:\20140320022038047.jpg");
            Emgu.CV.CvInvoke.cvSmooth(image, image2, SMOOTH_TYPE.CV_GAUSSIAN, 5, 5, 9, 9);
            Emgu.CV.CvInvoke.cvAddWeighted(image, 1.5, image2, -0.5, 0, image);

I know this answer is very late but I think this question worth an answer

The algorithm you used goes as follows:

Blur the original image using the Gaussian filter with given Mask Size & Sigma
Subtract the blurred image from the original (result is called Mask) to eliminate background and get the edges regions
Add a weighted portion from the mask to the original image by multiplying the Mask (the edges only) by K to enhance edges regions

K: User input (If K = 1 Unsharp, If K > 1 Highboost)

Mask(x,y) = Orig(x,y) – Blurred(x,y)

Result(x,y) = Orig(x,y) + K × Mask(x,y)

Here is the code:

public static Image<Gray, byte> Sharpen(Image<Gray, byte> image, int w, int h, double sigma1, double sigma2, int k)
    {
        w = (w % 2 == 0) ? w - 1 : w;
        h = (h % 2 == 0) ? h - 1 : h;
        //apply gaussian smoothing using w, h and sigma
        var gaussianSmooth = image.SmoothGaussian(w, h, sigma1, sigma2);
        //obtain the mask by subtracting the gaussian smoothed image from the original one
        var mask = image - gaussianSmooth;
        //add a weighted value k to the obtained mask
        mask *= k;
        //sum with the original image
        image += mask;
        return image;
    }

I think your problem is that you didn't subtract the smoothed image from the original image, anyway I hope this help


MATLAB: get members of histogram populations EXCLUDING those inbetween The Problem

I have a (noisy) square wave signal and need to extract the data points that make up the 'baseline' and 'plateau' ONLY. I do not want data from the rise and fall. So far I have accomplished the following:

data points designated as 'baseline' or 'plateau' superimposed onto original signal

I have defined 'baseline' data points as those from the minimum to the mean of the lowest 25% of signals + 1 standard deviation of the lowest 25% of signals.

Similarly I have defined 'plateau' data points as those from the maximum to the mean of the highest 25% of signals - 1 standard deviation of the highest 25% of signals.

This works ok, but as you can see: sometimes 1 standard deviation is not enough to encapsulate all of the baseline noise. 2 Standard deviations often encapsulate too much of the rise/fall of the signal.

Can anyone provide a more elegant solution to this?

My aim is to measure the signal to noise ratio of the square pulse, so I need to extract the populations that will give me baseline noise and the mean pulse height.

Assumption

Your noise is Gaussian.

Proposed solution

Compute the histogram of your data. There will be two peaks, as in the histogram you've shown. Extract the approximate positions of the peaks and use these values as starting values to fit two Gaussian curves to both peaks, each with the fit parameters $\mu$ (mean signal, i.e. position of the Gaussian) and $\sigma$ (width of the Gaussian).

Your signal will then be the difference in positions of the Gaussians, and your noise is given by their widths.

Yes, there will still be contributions from the sides of the rect-pulse. However, there probably is no way to fully get rid of them, since you need some criterion to tell apart plateau values from values from the side lobes, what always becomes problematic if a point is close to the baseline/plateau. You could mitigate this effect by using the histogram values not only as datapoints for the fit, but also as fit weights. By doing so, the (low number of) side lobe values gets less weight compared to the (high number of) plateau/baseline values.


Baseline Correction: What is the concept of a baseline shift and baseline correction? I have a laser induced fluorescence signal data-set. Currently the processing pipeline performs baseline correction by manually moving the signal so that the initial point is at the origin. Is this a right approach?

Is it to be done only if the shift is predetermined or known?
Does that mean baseline error signal has zero frequency?
Or can there be low frequency components which can be called baseline errors?

Edit:
Laser Induced Fluorescence signal is used to get relative quantities of collagen, NADH etc from the tissues.

Raw signal:

I am from neither signal processing nor physics background. The requirement I have been given is to develop a software for analysing the preprocessed data. However, I observed that currently the preprocessing uses the above said manual method for baseline correction.

Baselines are visually simple, most of the time, but can be a pain to filter.
With your graph, for instance, what should happen between 380 and 400 is not clear.

With fluorescence, it is not sure you can get a fine model of the baseline. So effectively, in many analytical instruments, there are more or less automatic baseline fits: constant, linear, parabolic.

Without hard models, it is common to use soft ones, but we need to know what you want to quantify in your data. If it looks like other analytical chemistry data, where you want to measure peaks above a baseline, you may try the technique given in Chromatography Baseline Placement, combining low-pass filtering and peak preservation.

Is there any chance you could share the dataset?


Modulation classification via demodulation confidence? I'm not sure if this is the right place to post this question, but will try anyway :-)

I recently started reading up on the field of automated modulation classification (AMC) where an algorithm tried to classify which modulation scheme an RF signal is using (eg BPSK, QPSK, etc..). Some if the common approaches include higher order statistics, deep neural networks, etc...

My question is, at least for digital modulation schemes with constellation diagrams, why cant this problem be solved by a very simple solution: try decoding the signal using various modulation schemes and see which gives you the most confidence. Confidence can be measured by something like mean distance to the closest constellation point of that scheme.

I realize that this may not be the fastest solution to the problem, and it cannot easily be applied to analog schemes (e.g. FM), but would it at least work? What would be its biggest issues?

Thanks

(PS: my background in not anywhere close to digital communications, so forgive any ignorance)

try decoding the signal using various modulation schemes and see which gives you the most confidence

well, simply because your "confidence" measure is basically equivalently hard to design generally enough as the other methods you mention.

Also, the search space is simply too humongous. How would you decode all possible combinations of constellations, symbol rates, with unknown phase and timing offsets, over all possible combinations with multi-carrier systems without any prior attempt to narrow things down? How low would your detection thresholds have to be so you could classify even without having done system-specific synchronization (which you can't do, without first assuming a system)?


Frequency-specific Delay? So an equalizer manipulates the gain of certain frequencies by use of a variety of filters, as we all know. But is there a way to manipulate delay in a frequency-specific manner?

A parametric equalizer

Above is a link to a picture of a parametric equalizer. What if the white line in the picture represented delay time instead of gain?

Is this possible? Does it exist? And if so, what is it called?

Yes it does exist and these are called "all-pass" filters in that over the band of interest the magnitude does not change (other than a possible fixed gain at all frequencies) but modify the phase response. The time delay at a specific frequency is the negative derivative of phase with respect to frequency, so the goal of these filters is to selectively change the slope of the phase according to the desired delay at each frequency.

The primary characteristic of a digital all-pass filter is that it will only have poles inside the unit circle on the z-plane and for each of those poles will have a conjugate reciprocal zero outside the unit circle, as depicted in the plot below.


Characteristic and moment generating function of a random variable interpretation I have been studying about moments and cumulants of a random variable. Even though the definitions of characteristic and moments generating function are very similar (only the sign in the exponential is changed) to Fourier and Laplace transforms, respectively.

$\Phi_x(\xi) = E [ e^{j \xi x(\zeta)} ] = \int_{-\infty}^{\infty} f_x(x) e^{j \xi x}  dx $

$\Phi_x(s) = E [ e^{sx(\zeta)} ] = \int_{-\infty}^{\infty} f_x(x) e^{sx}  dx $

The book I've been reading says the variable ξ is not and should not be interpreted as frequency!

How is it supposed to be interpreted? The functions $\Phi_x(\xi)$ and $\Phi_x(s)$ possess all the properties associated with the Fourier
and Laplace transforms, respectively.

Were these exponentials placed there just to be able to perform the Taylor series decomposition and then take the mth-order derivative of every term and extract the appropriate mth-order moment ?

As you well said, these functions are useful just for their mathematical properties that allow to compute the m-th order moment of a random variable.

Additionally, they happen to match Fourier transform theorems. However, not in the sense of a time-domain signal, since the variable in the probability density function (pdf) is the random variable itself, and not time.
A frequency interpretation only makes sense in time-dependent functions. Actually the definition itself of frequency is how often something happens.

Anyway, one could consider the "bandwidth" for $\Phi_X$ in the $\xi$-domain as how spread the pdf is around the mean.
A more spread pdf would imply a narrower "bandwidth", and less spread of the pdf would imply a wider "bandwidth".
For example, consider a uniformly distributed random variable:
$$f_X(x) = 1/\delta~\Pi(x/\delta)$$
$$\Phi_X(\xi) = sinc(\delta~\xi)$$
As $\delta$ gets smaller, the random variable is constrained in a narrower interval, while the "bandwidth" in $\Phi_X(\xi)$ gets wider, and vice versa (measuring the "bandwidth" as the width of the main lobe, $2/\delta$, as we do for communication signals).

Just by observing the characteristic function one may have an idea of how spread is the pdf. By comparing two characteristic functions with different "bandwidths" one may determine which pdf is spread the most.


New Horizons at Pluto I'm teaching an undergraduate linear algebra course, and looking for applications.  The New Horizons spacecraft is approaching Pluto, and returning images I find fascinating:
http://pluto.jhuapl.edu/index.php

They briefly refer to 'image deconvolution'.  I expect they are using the most sophisticated tools available, some nonlinear, but my vague understanding is that some simple methods are linear algebra based.

Can anyone provide a reference or expository paper on linear algebra methods for image deconvolution that would be a good starting point?

When talking about image deconvolution, it is referring to the process of correcting an image for distortion caused by the optical properties of the camera or imaging system. This distortion is characterized by the point spread function of the optical system. The point spread function describes the response of an imaging system to a point source. For a non point source ,the image is a superposition of the point spread function of multiple point sources representing the non point source.

There are many techniques of image deconvolution and references, too many to list in this format. But for examples of applications using linear algebra, you may want to begin by looking for references on image deconvolution using
Wiener filters, or the Richardson-Lucy method. Both require knowledge of the point spread function of  the system , but both involve linear algebra.


Convolution - Massive spike in noise section I am computing the correlation with MATLAB using the convolution fonction with a signal flipped before feeded into the function. Here are the results.

The first plot is the actual signal, the second one is the signal I am trying to find in the first signal and the third one is the convolution results. On the first plot, we can see at $n=3000$ and more  that the signal we are searching for is there and it is confirm by the spike in the third plot at $n = 3370$. However, there is massive spike in the beginning of the convolution (plot #3/ $n < 3000$) even if the beginning of the first signal is noise (plot #1/ $n < 3000$) . Why is there such spike in the convolution plot.

If your searched-for signal is partially similar to itself at some offset, that can also show up when correlated against an actual input signal containing it.

Test your searched-for signal by autocorrelation against itself to check for this problem.  If you want a sharp convolution peak against test input, you need a searched-for signal that has very low self-similarity at any lag.


Distortion of OFDM signal when using RRC filter I am trying to transmit OFDM signal in pass-band based on RRC filter, but what I noticed that the signal before the RRC filter is completely different compared with signal after the RRC filter (It means the signal is destroyed when using RRc filter).  I tried to follow the signal without using any channel, I noticed the issue starts when using the RRC filter. For example,  I used only the system as follows:
signal --> ifft --> upsampling --> RRC filter --> RRC filter --> downsampling --> fft --> received signal

What is expected is the received signal equals to signal, but I get very different results.
what I get is a different signal. Why ? should that be equal? I think I must be the same.
Finally , when I delete the filter, Everything is OK.

I checked your code in my PC, you need just to delete the delay added before the filter. For example, you can use:
U_aft_fil = U_aft_fil(fil_delay+1:end);

Then when filtering it again at the receiving side, you delete it again :
U_r_fil = U_r_fil(fil_delay+1:end);

Good luck


Cosine of amplitude 1 but GNU Radio shows FFT has amplitude > 30,000 I am doing an FFT on a cosine wave which has an amplitude of 1. The FFT is amplitude is not 1 but over 30,000. What am I missing here?
It's from the GNU Radio FFT wiki page:

Consider the formula for the DFT (which the FFT efficiently computes as an algorithm):
$$X(k) = \sum_{n=0}^{N-1}x(n)e^{-j2\pi nk/N}$$
Notice that it is a summation over $N$ samples total. Also note using Euler's formula that a cosine function can be expressed as two exponential terms (that are visible in the DFT result) as:
$$x[n] = cos(2\pi f n + \theta) = \frac{1}{2}e^{j(2\pi f n + \theta)} + \frac{1}{2}e^{-j(2\pi f n + \theta)}$$
Observe how in the DFT equation that we multiply the two exponential terms in x[n] by $e^{-j2\pi nk/N}$, that latter which is a phase rotation (the exponents which is the phase add in the product, and generally the term $e^{j\phi}$ is simply a complex vector with magnitude $1$ and angle $\phi$ just in case there was any confusion in that; $e^{j\phi}$ is simply a phase rotation). For each of those two exponential terms, there is a value of $k$ that would cause one of them to rotate to the same phase for every $n$, which then would grow to the large value in the summation. Given the OP's summation has grown to approximately 34,000, and that the magnitude of each of the exponential terms is $1/2$, suggests that there are approximately 68,000 total samples in the DFT assuming no windowing is also applied to the signal (which would then require even more samples due to the amplitude reduction through windowing). For every bin in the DFT result, this is the correlation to $e^{j\omega_k n}$ (a spinning phasor at frequency $\omega_k$, with $\omega_k = 2\pi k/N$), and sinusoidal functions have two such terms as given my Euler's formula: The worlds should now align.
$f$ is the fractional frequency term typically given as a value between $0$ and $0.5$ for the form in the equation above, where it represents the continuous time frequency divided by the sampling rate. Here we see further the values for k where that occurs:
When $k = Nf$
$$X(k=Nf) =  \sum_{n=0}^{N-1}x(n)e^{-j2\pi fn}$$
$$ = \sum_{n=0}^{N-1}\bigg( \frac{1}{2}e^{j(2\pi f n + \theta)} + \frac{1}{2}e^{-j(2\pi f n + \theta)} \bigg)e^{-j2\pi fn}$$
$$ = \frac{1}{2}\sum_{n=0}^{N-1}e^{j\theta} + \frac{1}{2}\sum_{n=0}^{N-1}e^{-j(4\pi fn + \theta)}$$
The first term grows to $\frac{N}{2}e^{j\theta}$ in the summation while most of the individual phasors in the second term will cancel leaving either zero or small residual depending on the value for $f$.
Similarly consider when $k=N(1-f)$
$$X(k=N(1-f)) =  \sum_{n=0}^{N-1}x(n)e^{-j2\pi(1-f)n}$$
$$ = \sum_{n=0}^{N-1}x(n)e^{-j2\pi n}e^{+j2\pi fn} = \sum_{n=0}^{N-1}x(n)e^{+j2\pi fn}$$
$$ = \sum_{n=0}^{N-1}\bigg( \frac{1}{2}e^{j(2\pi f n + \theta)} + \frac{1}{2}e^{-j(2\pi f n + \theta)} \bigg)e^{+j2\pi fn}  $$
$$ = \frac{1}{2}\sum_{n=0}^{N-1}e^{4\pi fn + \theta} + \frac{1}{2}\sum_{n=0}^{N-1}e^{-j\theta}$$


FFT frequency question Sorry if this is a stupid question but so far i have learned that after getting the results of a FFT (an array of complex numbers) you can get the frequency of a sound sample via the following:

1) Calculate the peak for each complex number: $\sqrt{\mathrm{real}^2 + \mathrm{imag}^2}$

2) Get the index of the array element with the highest peak.

3) Index $\times$ mhz of the sample audio $/$ number of elements in the array

This seems to work however I'm wondering if the second-highest peak would be the second "most present" frequency, and so on. So for example, if I had a guitar chord, and I got the top 6 frequencies, would those be the notes of a 6 stringed chord?

Unlikely.  The top 6 FFT magnitudes would more probably contain several overtones or harmonics that are stronger than the 6 note pitches in a chord.  Furthermore, some of the spectral frequency peaks can be between FFT bins, depending on the length of the FFT.  So the frequency you've calculated by just multiplying by a bin index could be different from any pitch frequency in the audio.  Depending on any window used, you might also find windowing artifacts in the FFT magnitude spectrum.


Confusion with regards to units in Sinc pulse I am currently reading a section in the Handbook for MRI Pulse Sequence book where it says that the time-bandwidth product for a SINC pulse is given by:

$$
T \Delta f = Z
$$

where $T$ is the length of the pulse, $\Delta f$ is the bandwidth and $Z$ is the number of zero crossings of the SINC pulse. Now, my question might be very naive but it says that this is a dimensionless quantity but I do not understand why that will be:

For example, if the bandwidth is described in radians/sec, does it not have units of radians? In that case, how does it relate to the number of zero crossings, which is indeed a discrete number.

The length of the pulse, $T$, would be in units of time (for example, in seconds). The bandwidth would be in units of frequency (for example, in Hz). So

$$[T][\Delta f]=sec\frac{1}{sec}=1$$

That means that $Z$ is a dimensionless quantity. That makes sense, as $Z$ is just a natural number that describes the number of times that something happens (in this case, the number of zero crossings), and should not have units at all.


Discrete-time, second-order section, direct-form II filter I found dfilt.df2sos in Matlab and trying to figure it out how do I do that in Python? Is there any function that works similarly with the Matlab function?
I also tried to browse but I could not find a similar function in Scipy.

More dfilt.df2sos https://www.mathworks.com/help/signal/ref/dfilt.df2sos.html

Best,

This is Direct Form 1 but it's still cascaded SOS.  it's in C, but you should be able to translate it to Python easily.  the register keyword is not important.
typedef struct
{
float filterCoefficients[5];    /* b0, b1, b2, a1, a2   (a0 is normalized to 1) */
} sectionData;

void filterBuffer(int32 Nsamples, int Nsections, float* y, float* x, float* filterStates, sectionData** sections)
{
sectionData** section_ptr;
float sampleValue, state1, state2, *state_ptr, *coef_ptr;

for (int n=0; n<Nsamples; n++)
    {
    section_ptr = sections;                     // reset section pointer
    state_ptr = filterStates;                   // reset state pointer

    sampleValue = x[n];                         // get input sample

    state2 = *state_ptr++;                      //   x[n-2]
    state1 = *state_ptr--;                      //   x[n-1]
    *state_ptr++ = state1;                      //   x[n-1] -> x[n-2]
    *state_ptr++ = sampleValue;                 //   x[n]   -> x[n-1]

    int i = Nsections;
    while (--i >= 0)
        {
        coef_ptr = &((*section_ptr++)->filterCoefficients[0]);  // point to section filter coefficients

        sampleValue =  *coef_ptr++ * sampleValue;   //   b0*x[n]
        sampleValue += *coef_ptr++ * state1;        //   b1*x[n-1]
        sampleValue += *coef_ptr++ * state2;        //   b2*x[n-2]

        state2 = *state_ptr++;                      //   y[n-2]
        state1 = *state_ptr--;                      //   y[n-1]
        *state_ptr++ = state1;                      //   y[n-1] -> y[n-2]

        sampleValue += *coef_ptr++ * state1;        //  -a1*y[n-1]
        sampleValue += *coef_ptr++ * state2;        //  -a2*y[n-2]

        *state_ptr++ = sampleValue;                 //   y[n]   -> y[n-1]

        //     at this point sampleValue is the output of the section just completed
        //     and will be the input to the section about to be computed unless this
        //     is the last section, in which case sampleValue is the output.
        //     state1 and state2 are the output states of the section just completed
        //     and will be the input states of of the section about to be computed.
        }

    y[n] = sampleValue;                             // store output sample
    }
}


Integral of an upsampled signal, without actually resampling it Given a signal X which is sampled at a certain frequency. The value we currently compute is given as the integral of the upsampled signal. Thus: Y = X but 100 times upsampled, by means of sinc interpolation or by using an FFT resampler. The integral is simply the sum of all values in Y.

The calculation of this integral is easy, yet I would like to speed it up by avoiding the upsampling step. Is there any possibility to obtain the integral of the upsampled signal, without actually resampling it ?

Preserving integral through downsampling and Do lowpass filters affect the integral over the signal? are similar questions, but deal with cases of downsampling or bandlimiting. In both cases, it is clear that an integral will not exactly be preserved by downscaling.

Yet, in this question, we are talking about upsampling. I suspect this is possibly because X does contain all information necessary to create Y.

If you upsample a causal discrete-time signal $x[n]$ by an integer factor $L$ and you use an interpolation filter with a causal impulse response $h[n]$, the upsampled and interpolated signal is given by

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

(Note that the requirement of causality can be easily removed; I just use it to have the lower summation indices start with the value $0$.)

If you want to compute the sum of $y[n]$, you would need to compute

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

However, you can interchange the sums in (2), resulting in

$$s_y=\sum_{k=0}^{\infty}x[k]\sum_{n=0}^{\infty}h[n-kL]\tag{3}$$

The value of the inner sum in (3) is independent of the index $k$, so we can write

$$s_y=s_h\cdot \sum_{k=0}^{\infty}x[k]=s_h\cdot s_x\tag{4}$$

with

$$s_h=\sum_{n=0}^{\infty}h[n]$$

So the sum of the upsampled and interpolated signal is simply the sum of the original signal times a factor which is given by the sum over the impulse response values of the interpolation filter. Note that for the ideal sinc filter, this factor equals the upsampling factor $L$. This is also at least approximately true for all practical interpolation filters, because the sum over the impulse response equals the DC gain of the interpolation filter, which should ideally equal $L$.

So if you want to keep things simple (and why wouldn't you?), you could use

$$s_y\approx L\cdot s_x$$

which is a very good approximation (or even exactly true) for any reasonable interpolation filter.

In the above derivations I've assumed that all the (infinite) sums exist, which in practice is always the case if we deal with finite length signals and with stable interpolation filters.


FFT phase spectrum to suppress noise? All that I heard so far about phase spectrum - that it is very hard to process. Are there any methods to do that, and is it possible to extract information about the noise from phase spectrum?

you can look at

https://math.stackexchange.com/questions/1290620/what-is-the-importance-of-phase-spectrum-in-fourier-transform

the short answer is that for noise, phase is random.  There shouldn't be any structure to exploit, while for signal, there often is.


Texture mapping on a 3D face from 2D face image I have 3D points (xyz) of a face image and I have 2D face image of the same person.

The 3D points of the face are such that If I project the points to 2D plane, it matches with the 2D face image.

I want to map the color information on to the 3D image.
So after mapping I will have the texture of each point.
I want to do it using Matlab/OpenCV.

Can anyone suggest any algorithm or refer some source for me to start.

What you're trying to do is more in the image synthesis / computer graphics domain (technically, it's "texture mapping") than it is in the image-processing domain, so it's not obvious that Matlab/OpenCV are going to have any relevant supporting functionality.

If I had the same problem in my (preferred) Python/SciPy/Matplotlib world, I'd probably see if mplot3d's surface plotting was up to the job (passing the 2D image as the facecolors parameter).  However general experience is that it doesn't deal well with 100s x 100s size grids, only a few 10s x 10s.  Sorry, I don't know what Matlab's equivalent is.  (And if mplot3d wasn't up to the job I'd be looking at the mlab VTK wrapping's texture mapping support).


Special case: Band-limited in time domain and band-limited in frequency domain I know except for some special cases, aliasing is unavoidable. Assume we time-limit a function, $f(t)$, so that it is zero outside an interval say $[0,T]$ to form $y(t)$. Then, in the frequency domain, infinite frequency components are introduced.

If the system that makes this time limitation has impulse response $h(t)$, we have its Fourier representation $H(u)$. We can then calculate $Y(u)$ using convolution to obtain:
$$Y(u)=F(u)*H(u)$$
Obviously, it has frequency components that are replicated to infinity. This phenomenon is because of sliding it across the original function to calculate Convolution.

Could anyone possibly let me know in what condition we can have band-limited input and output.

The condition is that the function $f(t)=0$, if you are talking about continuous time signals and frequencies (some of your formulations are confusing to me). In other words, a (nonzero) time-limited signal cannot be also band-limited. In other words, a function and its (continuous) Fourier transform cannot both have finite support.

The core of one proof (using complex analysis) relies on Paley-Wiener theorems: the Fourier transform of a time-limited function would be an entire function (a complex-valued function, holomorphic over the whole complex plane), that would vanish on an open interval, and hence be zero (almost) everywhere.

You can find details in Fourier transforms of compactly supported functions, and the paper (for instance) The uncertainty principle for Fourier transforms on the real line. This is a reason for the development of time-frequency analysis.


Variance of correlated noise I think this is an easy question but I am confused about the answer. Assume we have correlated noise distributed as $$n \sim \mathcal{N}(0,\sigma^2 \:A),$$ where $A$ is a given square matrix. If we multiply this noise by $A^{-1}$, i.e., $$w = A^{-1} n.$$ Then, the noise covariance matrix in $w$ is given as $$\mathbb{Var}\{w\} = \mathbb{Var}\{A^{-1} n\} = (A^{-1})^2 \mathbb{Var}\{n\} = (A^{-1})^2  A = A^{-1}.$$

Am I missing something?

Even though your calculation yields the correct result (in this case), the steps are not completely correct.

First of all, the covariance matrix of a random variable $n$ is given by

$$ C = E[nn^H]. $$

The nomenclature of variance is more suitable for scalar random variables, in my opinion (or it refers to the diagonal of the covariance matrix, i.e. the variance of each element of $n$).

Now, you state the covariance of $n$ should be $\sigma^2A$. Hence, we have

$$ E[nn^H] = \sigma^2A $$.

Now, you are interested in the covariance matrix of the random variable $y=A^{-1}n$, which is given by

$$
D = E[A^{-1}nn^HA^{-H}]
$$
Now, knowing that $A$ is constant, we can take it out of the expectation, to get

$$
D = A^{-1}E[nn^H]A^{-H}=A^{-1}\sigma^2A A^{-H}=\sigma^2A^{-H}.
$$

Now, knowing that $A$ is a covariance matrix and these are generally semi-positive definite and symmetric, we have $A^{-1}=A^{-H}$, so the final result is

$$
D = \sigma^2 A^{-1}
$$

which matches your guess. Note that in case you multiply $n$ with a different matrix $B$, the result is not $\sigma^2B^2A$ (as you might have guessed), but according to above calculation it is $\sigma^2BAB^H$.

I often refer to the following link, when it comes to expectations that involve matrices: http://www.psi.toronto.edu/matrix/expect.html


How would I calculate the DC Gain for a Chebyshev Type 2 filter? I'm currently working on an implementation of a Chebyshev Type 2 filter but the gain is all over the place when playing around with the order and ripple. I haven't found anything online but (as far as I understand) I need to calculate the DC gain to compensate for gain differences. How would I go about doing that?

I've tried to do s = 0 for my analog prototype which leaves me with just the coefficients b1 / a2 for every biquad. I multiplied all the fractions of the individual biquads together but the result did not help at all, it seems just as random as the gain itself.

Can anyone point me in the right direction?

I suppose that you obtain a transfer function of the desired discrete-time system in the form

$$H(z)=\frac{b_0+b_1z^{-1}+\ldots +b_Nz^{-N}}{1+a_1z^{-1}+\ldots +a_Nz^{-N}}\tag{1}$$

From that transfer function you compute the coefficients of the second-order sections. Before doing that, you can make sure that the DC gain of $(1)$ equals $1$ (which of course only makes sense for a low pass or a band stop filter).

The DC gain of $(1)$ simply equals

$$H(1)=\frac{b_0+b_1+\ldots +b_N}{1+a_1+\dots +a_N}$$

Obviously, the DC gain is independent of the chosen filter type.


Why we add imaginary part in inverse fourier transform since the time domain signal has only real value? The equation of inverse Fourier transform is the following:
$$ f(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} F(\omega) e^{i\omega t} d\omega $$
$$ f(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} F(\omega) \big[\cos(\omega t) + i\sin(\omega t)\big] d\omega $$
Now, if I have only real value in the time domain $f(t)$, then the imaginary part of $\int e^{i\omega t}$ must go to zero. In our signal processing course, I always had real valued time domain function. Then why am I adding this imaginary part of $\int e^{i\omega t}$ in inverse Fourier transform since it always goes to zero?

then the imaginary part must go to zero

No. That's a fallacy. There's really not much to explain here other than this intuition is wrong.
The imaginary part of a Fourier transform is only zero if the transformed signal has conjugate symmetry (which is the same as symmetry for real signals) around t=0.
Think about it. You decomposed your $e^{-ix}=\cos x + i\sin x$, and then used an integral to project your signal onto these two things. (since the integral of a product between a signal and another (conjugated) signal is the inner product of the signal vector space, and hence the projection, if you're after the mathematical formalism.)
The part of your signal that is symmetrical to t=0 projects onto the cosine, and the part that is antisymmetrical to t=0 to the sine. The sine has an imaginary unit in front.
This is really among the first things that are typically illustrated in textbooks when they introduce the Fourier transform, so this might be a great point to go back to that!


Converting IQ to RF data I appreciate any help: I wanna convert IQ data to its RF data, representing an image. The center frequency and sampling rate are available. What I do is:

First I do up-sampling in MATLAB by applying interpolation along depth. So, the new data has a new sampling frequency higher than the first one.

My questions are: 1) do I need to do low-pass filtering before multiplying the data by the exponential term (that converts IQ to RF)? 2) what is the exponential form that I need to use to multiply the IQ data by?

I was wondering if the length of one of my vectors is, say, 1000, ( Calling it vector A), then how to do zero padding? (I assume a number of zero elements are added to both start and end of this vector and then the Low-pass filtering is done) How to decide about the number of zeros added? Sorry, if these seem to be too simple questions

I've read the documentation I copy here and it is confusing for me: (Page 8 and afterwards)

http://folk.ntnu.no/htorp/Undervisning/TTK10/IQdemodulation.pdf

Thanks

For ideal performance, you can DFT your input, pad out the spectrum with zeros, and inverse DFT the result.  This is not a particularly efficient approach, though.

Pretty much the same result is obtained with interpolation by zero-padding and then running that stream through a low pass filter.  An FIR sinc filter is preferred for this for two reasons--(1) it has nearly ideal performance and introduces negligible distortion into the baseband signal, and (2) when zero padding, most of your input samples are zero.  Since 0*x=0, you don't have to do all those multiplies; for your 6:1 example, you can design a 72-point FIR and still only do 12 multiplies per output sample.

Once the filtering is done, you can upconvert (multiply by the exponential).  This should be $e^{j\omega t}$ where $\omega=2\pi f$ and f is your RF carrier frequency.


Generalized Likelihood Ratio test for correlated data Consider a sequence of random variables $\overline{z} = \{z(k-M+1), ..., z(k)\}$, with probability density depending upon a scalar parameter $\theta$. It is intended to decide between two hypotheses

$$
H_0: \theta = \theta_0    \hspace{3em} for \hspace{3em}  k-M+1 \leq i \leq k \\
H_1: \theta = \theta_1    \hspace{3em} for \hspace{3em}  k-M+1 \leq i \leq k \tag{1} $$

It is assumed that all the parameters are known except for $\theta_1$ which is replaced by its MLE estimate.
Then, the GLR test is

$$L(\overline{z}) = \frac{p\left(\overline{z}, \; \hat{\theta_1}; \; H_1\right)} {p \left(\overline{z}, \; {\theta_0}; \; H_0 \right)} \tag{2}$$

where $\overline{z} = \{z(k-M+1), ..., z(k)\}$ is iid.

Then GLR test can be simplified to

$$L(\overline{z}) = \max\limits_{\theta_1} \left(  \sum\limits_{i = k-M+1}^{k} ln \frac{p(z(i)|\theta_1)}{p(z(i)|\theta_0)}  \right)\tag{3} $$

But what if $\overline{z} = \{z(k-M+1), ..., z(k)\}$ are correlated. I need to simplify $(2)$ for simplified Gaussian process.

Any help, referred references for correlated GLR test is highly appreciated.

A simple approach would be to transform your observations $\bar{z}$ into independent random variables is by using whitening transformation, i.e., if the covariance matrix $\Sigma$ is known, then find $W$ such that $W^TW = \Sigma^{-1}$. Next, define a new random variable $\bar{y}$ as $W \bar{z}$. The above transformation is only possible when $\Sigma$ is invertible.


How to downsample a matrix in columns ? MATLAB I am studying wavelets in image processing. I would like to learn how to downsample a matrix in columns using MATLAB. I have used the downsample(x,n) command to downsample the given matrix in rows. But how to do the same in columns?

x =

 1     2     3
 4     5     6
 7     8     9
10    11    12

y = downsample(x,2)

y =

 1     2     3
 7     8     9

But I would like to get the following result.

y=
     1       3
     4       6
     7       9
    10      12

How to do it in MATLAB? Please help.

Maybe Off-Topic here, but anyway !

Try it

y = downsample(x',2)'


Spectral analysis of a cross correlation function So it's widely known that the Cross Correlation of 2 signals helps us in figuring out the time delay in those signals by analyzing the peak of the correlation coefficient in the time domain.

For something I am working on, due to a short baseline in comparison to the wavelength of the signal in between two of my receivers, the cross correlation of the signals from my two sensors give a peak at zero which makes it hard for me to analyse them.

Which is what brings me to my original question i.e
What if I analyse the fft of this correlation coefficient array? Will I find anything hinting towards the time delay?

If not, is there anyway to figure out how I can find the time delay with such a short baseline (~50m) for a long wavelength (~20km) signal?
Thanks for the help. Any little help would be appreciated.

cross-correlation would determine the time-delay. but the scale of signals' length and the delay value should be considered precisely. I mean if you have two long signal of, for say, 100sec and you are looking for a time delay in an order of 1msec, you don't need to do cross-correlation over the whole 100sec. Instead, choose a proper window time (in msec order) of two signal and then perform the cross-correlation.


Gaussian White Noise - Relation Between Distribution and Correlation Im a beginner in signal processing so my question may be obvious.

A white noise has the property to have its autocorrelation function that is equal to

$$\mathbb{E}[f(t+\tau)f(t)]=\sigma^2 \delta(\tau).$$

We say that random signals can be Gaussian white noise, but for me if the signal has a Gaussian distribution it is necesseraly a white noise because as soon as we know the distribution then $\mathbb{E}[f(t+\tau)f(t)]$ is automatically determined.

Am I wrong ?

When we say the stochastic process $X_{t_1}X_{t_2}\cdots X_{t_i}$ is a Gaussian Process with mean $\mu$ and variance $\sigma^2$ it only means
$$X_{t_i}\sim\mathcal{N}(\mu,\sigma^2)$$
So still we know nothing about the relation of the $X_{t_i}$'s with each other.
When we add that the process is also white, it means the $X_{t_i}$'s are uncorrelated. That is

$$\mathsf{E}[(X_{t_i}-\mathsf{E}X_{t_i})\mathsf{E}(X_{t_j}-\mathsf{E}X_{t_j})^*]=\begin{cases}\sigma^2,&i=j\\0,&i\neq j\end{cases}$$

Off topic, but notice that the $\tau$ notation is an implicit assumption of the process being stationary.

The process can be Gaussian but colored (rather than white). In such case, the random variables are Gaussian but correlated. For instance, the autocorrelation function in case of stationary Gaussian process you mentioned can be something like this:
$$R_{XX}(\tau)=\sigma^2e^{-|\tau|}$$


Interpolated FIR filter group delay I'm trying to design a digital low pass filter with a narrow transition band. My sampling rate is 25 kHz, the cut off frequency is 60 Hz & the transition band width is 4 Hz. I'm looking for about 40 dB attenuation in the stop band & 0.1 dB in the pass band.

I looked at the traditional FIR methodology using the Windowed Sinc filter but the number of taps required to get such a transition band is too high. Going by the reference here, we get a value of 25000 which is 1 second in terms of time.

As an alternative, I looked into IFIR filters. While the model filter results in a decrease in the filter length by a value that's equal to the interpolation factor, the interpolated filter formed by inserting zeros into the model filter's kernel increases the length to almost the same as what you'd get with a conventional FIR filter.

So, my question is, does the IFIR methodology make the filter more efficient from a group delay point of view? Is the only benefit here the fact that more zeros implies fewer multiplications & additions? I'm trying to achieve around 1 to 2 cycles of delay (in terms of time for a 50 Hz signal). Considering that the sampling rate is so high, shouldn't it be easier to achieve this? If not, what other methods could be used to obtain such a filter?

If you need to discriminate between 60 and 65 Hz signals, the fact that you're sampling at 25 kHz means that you're basically doing real time signal processing: the Nyquist Frequency is far above the frequencies of interest.

Thus, from classical Fourier theory, you're going to need on the order of 1.0/5 Hz = 0.2 seconds of data, no matter how you process it.  If you do FIR, that means that you'll need at least 5000 coefficients --- and probably more like the 25000 that you mention.

It's possible that you could use much smaller IIR filters.  You'll do less math, but you'll still have to chew through at least 5000 samples for fundamental Fourier reasons.  And you'll need to be careful with that IIR filter, because the pole will be very close to the unit circle: you'll probably want to convert your data to float or double rather than an integer format (which would likely work for FIR).

Now, if there is some prior knowledge you have about your signals --- "high SNR, knowledge of the phase of the 60 and 65 Hz signals..." then you can create discrimination algorithms that are faster.  But all you know is 60/65 Hz, you're sort of out of luck.


Plotting the frequency response of a FIR filter I am trying to figure out the mechanics of plotting the frequency response of a FIR filter. For example, I have used an algorithm (Parks-McClellan), to generate a low pass FIR filter with an even number of coefficients. I've got 34 tap values. Now I want to show the frequency response of this configuration in a normalized frequency 0 ... 1.0 (Fs) graph. The literature says, "The frequency response of the filter is computed by passing the array of coefficients through the discrete Fourier transform (DFT)." The text goes on to show this nice smooth graph, plotting magnitude vs frequency from 0.0 to 1.0. I would like to reproduce that graph.

So I sent it through my DFT code, well a 34 point DFT generates a very blocky DFT and it looks nothing like the graphic. I wanted more points of the DFT (between 0  and 34) and tried adding fractional points 0.0, 0.1, 0.2, .. 34.0) and got a very "loopy" DFT with lots of humps (and also not at all like the graphic), I tried "upsampling" adding zeros between the values, still no joy. I tried allowing the tap sequence to repeat, no joy. So I'm really curious how its done. Sending it to the MATLAB plot function produces the expected graph, so I have some confidence the data is good. Just my understanding of the mechanics of how to get the higher resolution graph are incomplete.

If you want to do it manually, it is not upsampled by adding zeroes between coefficients. It is usually done by adding extra zeroes after the real coefficients. When padded to some convenient power of two like 512 coefficients, FFT can be used. The resulting amplitudes need to be scaled due to the padding. Matlab and Octave have freqz function to plot this directly for you.


FFT algorithm gives wrong signs on my output numbers I was told from another stack exchange (stackoverflow and computer science stack exchange) to post here - so hopefully now i can get some help from people who understand the algorithm, but if this is off topic here then i am kinda running out of places to ask....
I wrote my own FFT algorithm to learn it - but i am getting a different result to what i was expecting (correct numbers but wrong signs), but i cannot figure out why.
I am hoping some one here who has a good understanding on the FFT algorithm can understand what i did wrong, i currently suspect its my twiddle factors - but i am not sure at this point...
Firstly i precompute all the twiddle factors and indices. The twiddles are calculated like this:
        for (int stage = 0; stage < _passes; stage++) //log(N)/log(2)
        {
            for (int i = 0; i < _N; i++) //N = 2^n
            {
                //wing span, 2^1, 2^2 etc
                float span = pow(2, stage+1);
                float k = (i * _N / span) % _N;
                float a = pi2 * k / _N;

                //iFFT uses complex positive coefficients: exp(i*2*pi*k/N)
                Vector2 twiddle = new Vector2(cos(a), sin(a));

                //FFT uses complex conjugate: exp(-i*2*pi*k/N)
                Vector2 twiddleConj = twiddle.ComplexConjugate();

                int arrIndex = stage * _N + i;
                _twiddlesR[arrIndex] = twiddle; //inverse array
                _twiddlesF[arrIndex] = twiddleConj; //forward array
            }
        }

This was the output for N=4:
  First stage:
  [0] + (1+0i)[2]
  [2] +(-1+0i)[0]
  [1] + (1+0i)[3]
  [3] +(-1+0i)[1]
  Second stage:
  [0] + (1+0i)[2]
  [1] +  (0-i)[3]
  [2] +(-1+0i)[0]
  [3] +  (0+i)[1]

I ran the FFT with this code:
for (int stage = 0; stage < _passes; stage++)
{
    for (int k = 0; k < _N; k++)
    {
        PerformFft(stage, k);
    }
}

My FFT function:
    private void PerformFft(int stage, int k)
    {
        int arrIndex = stage * _N + k;
        Vector2[] twiddleArr = _isForward ? _twiddlesF : _twiddlesR;
        Vector2Int twiddleIndex = _twiddleIndices[arrIndex];
        Vector2 twiddle = twiddleArr[arrIndex];

        //ping pong between two read/write arrays by checking if stage is even or odd
        Vector2[] writeArr = stage % 2 == 0 ? _complexArr2 :_complexArr1;
        Vector2[] readArr = stage % 2 == 0 ? _complexArr1 :_complexArr2;

        Vector2 complexA = readArr[twiddleIndex.x];
        Vector2 complexB = readArr[twiddleIndex.y];

        Vector2 result = complexA + complexB.ComplexMultiply(twiddle);
        writeArr[k] = result;

    }

For the output i get the correct numbers but i get the wrong signs, here is some data:
Input:
0 + 0i
1 + 0i
2 + 0i
3 + 0i
Forward FFT Output:
6 + 0i
2 - 2i
2 + 0i
2 + 2i

But according to this FFT calculator i am close to the correct answer but the signs are wrong: https://scistatcalc.blogspot.com/2013/12/fft-calculator.html
Can any one see what logic mistake i am doing here? Been stuck on this for a week so far.

Reversing the order of the input sequence would provide the same result as what the OP achieved, but not to say this is exactly why it is occurring here. Below shows the details of the DIT algorithm after each stage; comparing each element step by step with the code should reveal the actual error.

Interpreting the OP's results (I assume [x] represents the index) reveals the error is in the core 2 pt DFT butterfly operation.
While the 2 pt DFT is given as
$Y[0] = X[0] + X[1]$
$Y[1] = X[0] - X[1]$
While the OP is doing the following:
$Y[0] = X[0] + X[1]$
$Y[1] = X[1] - X[0]$
With $X[0]$ and $X[1]$ as the 2-Pt DFT input and $Y[0]$ and $Y[1]$ As the 2-Pt DFT output.
  First stage:
  [0] + (1+0i)[2]   = 0+2
  [2] +(-1+0i)[0]   = 2-0   (should be 0-2= -2!)
  [1] + (1+0i)[3]   = 1+3
  [3] +(-1+0i)[1]   = 3-1  (should be 1 -3= -2!)
  Second stage:
  [0] + (1+0i)[2]
  [1] +  (0-i)[3]
  [2] +(-1+0i)[0]
  [3] +  (0+i)[1]


A Better High-Order Low-Pass Filter? Does there exist a discrete low-pass filter algorithm that is greatly more computationally efficient for high-order (e.g. greater than 10) low-pass filters than straightforward methods like series of IIR Butterworth filters?

You can implement very efficient high order low pass filter using cumulative summing. The theory is as follows. A perfect integrator has the recursive definition

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

If you take the delayed difference of the output as in
$$y'[n] = y[n] - y[n-d]$$

you are realizing a boxcar filter of length $d-1$, as can be seen from explicitly resolving the recursion for $y[n]$ in $y'[n]$:

$$ y'[n] = (y[n-1]+x[n]) - y[n-d]$$

$$ y'[n] = ((y[n-2] + x[n-1]) + x[n]) - y[n-d]$$

$$ \dots$$

$$ y'[n] = y[n-d] + \sum_{k=1}^{d-1}x[n-k] - y[n-d] = \sum_{k=1}^{d-1}x[n-k]$$

That means you have a low pass filter of effective order $d-1$ from just one addition and one subtraction per sample.

Of course, boxcars are only useful up to a certain point. But the same idea can be generalized  to a larger class of filters with a trade off between filter shape and computational effort. And that is done by concatenating several integrator-difference stages.

Two in series will generate a triangular impulse response, three a parabolic one, and so on. Each stage convolves the output with yet another boxcar. In the frequency domain the original sinc() shaped response is taken to the n-th power (or if you have different delay lengths for the different stages, multiplied by other differently scaled sincs).

But both response sequences in time and frequency domain behave very predictably, because the central limit theorem of statistics predicts that the time domain response quickly approaches a Gaussian distribution. And because the fourier transform of a Gaussian is also Gaussian, you'll get a very good approximation of a Gaussian frequency response low pass filter after only a few stages.

One note about numerical precision though. Because you are integrating, the content of the delay buffer can be become very large and the difference precision suffers. This can be resolved by making the integrator leaky, so that it doesn't accumulate a constant function up to infinity. The leaky integrator difference equation is

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

for some sufficiently small positive $\epsilon$.


Issue with the time vector returned by $\tt signal.spectrogram$ function I'm puzzled by some very simple concept with building up a spectrogram.

Here is a toy example of the issue:

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

f0 = 5
t = 2
fs = 250
N = t*fs
t = np.linspace(0, t, N, endpoint=False)
x = np.sin(2*scipy.pi*f0*t)
f, t, Sxx = signal.spectrogram(x, fs)
plt.pcolormesh(t, f, Sxx)
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [sec]')
plt.show()

What I get for t puzzles me, I get [0.512,1.408]. I would expect to get an array starting from 0.0 to 2.0 seconds long - so if I want to see the power spectrum density in 1.8 seconds I could easily do so.

I probably miss something pretty basic, so it would be great if someone could shed some light on this.

Mike.

The default parameters of signal.spectrogram are:

nperseg = 256
noverlap = nperseg/8 = 32

This means that:

The length of analysis window is $256$ samples ($256/250 = 1.024$ second)
The overlap between consecutive windows is $32$ samples ($32/250 = 0.128$ second)

The timestamps returned by signal.spectrogram correspond to the centres of a window. So in your case, the first timestamp should be at $256/2 = 128$ samples, which is equal to $0.512$ seconds given a sampling frequency of $250 \mathrm{Hz}$.

The beginning of the next window will be shifted by $256-32 = 224$ samples ($32$ is an overlap), so its centre will be at $224 + 128 = 352$ samples (which corresponds to $1.408$ second). The end of the second window is at $480$ samples, so there is not enough data to process more frames.

This explains why your time vector is [0.512, 1.408].

Here is a graphical explanation:

So if you want to achieve a specific frequency spacing, you must modify the nperseg and noverlap parameters.


Fourier transform of a discrete periodic signal The given sequence
$$x[k]=\delta[k-1]-3 \delta[k] + \delta[k-2]$$
represents a period of some discrete periodic signal $e[k]$ with period N=3. Find the Fourier transform  of $e[k]$.

This is an exercise in my textbook and i don't get it. Is $x[k]$ the first period of $e[k]$? If so, does that mean that the DFT of $x[k]$ is the same as the DFT of $e[k]$?

For a discrete periodic signal you normally compute the Fourier series coefficients. They can be obtained from the DFT by dividing the output of the DFT by N, where number of samples in a period.

The DFT assumes the signal is N-sample periodic anyway.

In Matlab you will implement this as

x=[-3 1 1];
ck=fft(x)/3;


μ-law compression technique In μ-law compression formula i have always seen that the value of Xmax and Ymax is always same but can they be different?

Yes, they can be different, depending what your maximum input value is.


Fourier transform properties I have to find the Fourier transform of $$ x(t)= \frac{1}{T}e^{-\frac{t-T}{T}}u(t-T) $$

First I applied traslation property , so $$ F[x(t-T)] = X(f) e^{-i 2 \pi f T} $$

after I applied time scaling property , and I obtained
$$ F[x(\frac{t}{T})] = |T| X(Tf) e^{-i 2 \pi f T} $$

Now I calculate $$ \int_{T}^{+\infty } \frac{1}{T}Te^{-\frac{t-T}{T}}u(t-T) e^{-i 2 \pi f t } $$ that’s results $$ \frac{T}{1+ 2 i\pi f t } $$

So I obtained that the final Fourier transformation should be $$ X_f = |T| \frac{T}{1 + 2 i \pi f T} e^{- 2 i \pi f T } $$

On my book the result is the same except for T|T| but I don’t know what’s wrong. Thank you so much

If you already know the - quite well-known :) - FT pair $$e^{-at}u(t) \longleftrightarrow \frac{1}{a+j2\pi f}$$ then by setting $a = 1/T$, you get $$e^{-\frac{t}{T}}u(t) \longleftrightarrow \frac{1}{\frac{1}{T} + j2\pi f}$$
Time delay property $$x(t-t_0) \longleftrightarrow X(f)e^{-j2\pi ft_0}$$ can be used now by setting $t_0 = T$, and then $$x(t-T) = e^{-\frac{t-T}{T}}u(t-T) \longleftrightarrow \frac{1}{\frac{1}{T} + j2\pi f}e^{-j2\pi fT}$$
Adding your $1/T$ constant in front, you get $$\frac{1}{T}e^{-\frac{t-T}{T}}u(t-T) \longleftrightarrow \frac{\frac{1}{T}}{\frac{1}{T} + j2\pi f}e^{-j2\pi fT}$$
Finally, multiplying both numerator and denominator by $T$ you get what you want: $$\frac{1}{T}e^{-\frac{t-T}{T}}u(t-T) \longleftrightarrow \frac{1}{1 + j2\pi fT}e^{-j2\pi fT}$$


How Hessian feature detector works? I know about Harris corner detector, and I understand the basic idea of its second moment matrix,
$$M = \left[ \begin{array}{cc}
I_x^2 & I_xI_y \\
I_xI_y & I_y^2
\end{array} \right]$$,
 edges and other unstable points can be removed via $M$.

But about Hessian detector, it uses Hessian matrix to detect key points and remove edges,
$$\mathcal{H} = \left[ \begin{array}{cc}
I_{xx} & I_{xy} \\
I_{xy} & I_{yy}
\end{array} \right]$$,
and I don't understand how could $\mathcal{H}$ remove edge and detect stable points? What's the intuitive basic idea behind it?

I will try to avoid math, because math and "how to do it" tutorials can be easily found.

So, I start by pointing out one VERY important thing: One does not compute Harris for a single pixel, but for a vicinity (a patch of image) around that pixel! Let $I(i)_{xx}, I(i)_{xy} ...$ be your derivatives for a point $i_0$, then,

$H = \left[
\begin{array}{cc}
\sum_{i\in V}I(i)_{xx} w (i-i_0) & \sum_{i\in V}I(i)_{xy}w (i-i_0) \\
\sum_{i\in V}I(i)_{xy} w (i-i_0)& \sum_{i\in V}I(i)_{yy} w (i-i_0)\\
\end{array}
  \right] $

The $w(t)$ is a Gaussian kernel.
The previous eq tells you to integrate the derivative values over vicinity $V$ around current pixel. Each value of the neighbors is multiplied with a value that shrinks as the distance increases. The law of decreasing follows a Gaussian, because $w(t)$ is Gaussian centred at $i_0$. And that's it with math.

Now, back to the empirical observations. If you use solely the derivatives, and that pixel is part of a linear structure (edge), then, you get a strong response for the derivatives.
On the other hand, if the pixel is at a corner (an intersection of two edges) then, the derivative responses will cancel themselves off.

Saying that, the the Hessian is able to capture the local structure in that vicinity without "cancelling" effect. BUT very important, you have to integrate in order to get a proper Hessian.

Having a Hessian, obtained using Harris method or by other means, one might want to extract information about the vicinity. There are methods to get numerical values on how likely is to have an edge at current pixel, a corner, etc. Check the corner detection theory.

Now, about "stable points" or salient points. Picture that you are in a foreign town with no GPS and only with a good map. If you are "teleported" in a middle of a street, you might locate the street on the map, but you cannot tell where exactly are you on that street or in what direction you should go to move left or right (wrt to the map).
Imagine now that you are at a intersection. Then, you precisely can point your position on the map!. (Of course, assume that two streets doesn't intersect more than once).

Imagine now that you must match two images. One acts as a map, and the other as the city. You must find pixels that can be uniquely described, so you can do the matching. Check images on this post for example of matching. These points are called salient points. Moreover, the corner points tend not to change their 'cornerness' properties when the image is scaled, translated, rotated, skewed, etc. (affine transforms) This is why they are called "stable".

Some points in the image allows you to uniquely identify them. These pixels are located at corners or at intersection of lines. Imagine that your vicinity $V$ is on a line. Except for the orientation of the line, you cannot find anything else from that vicinity. But if $V$ is on a corner, than, you can find out the directions of the lines that intersect, maybe the angle, etc.

Not all corner points are salient, but only corner points have great chances of being salient.

Hope it helps!

p.s. How to find if a point is corner or not, take a look at Harris paper.

p.p.s. More on matching, search for SIFT or SURF.

p.p.p.s. There is a "generalization" of the Harris method, called Structure Tensor. Check Knutsson seminal work!


Distortion in phase difference measurement using FFT I am doing distance measurements by transmitting a 100 kHz tone and receiving it and measuring the phase difference.The distance is simulated by a delay simulator which is very precise (around 10 psec)
The 100 kHz tone is PM(Phase not PSK) modulated with an index of 1.1 rad and upconverted to 70 MHz. Then it is received downconverted and demodulated.
Everything seems fine except the phase difference varies within 0.15 degrees which is quite high as the tone has a high SNR >30 dB.
The steps leading to obtaining the phase difference are:

Sample the received 100 kHz tone at Fs= 10 MHz.
Take 32768 samples.
Apply Blackman-Harris window on the received samples.
Apply FFT with a length of 32768 on both the received samples and the reference tone.
Find the bin of the received and the reference tone: major_bin = floor(100000 * 32768/ Fs)=floor(327.68)=327
Find the phase of both the received and reference tones using:
atan(imaginary_fftsample[327]/real_fftsample[327])*180/pi
Subtract the phases to get the difference and calculate distance.

At this point I get variations in phase difference around 0.15 degrees while there is no change in the distance.
I suspect there is something wrong with the steps above. I don't want to use a PLL instead FFT.
Any ideas?

The OP mentions "demodulating" but I don't actually see that in the process (meaning the phase modulation is not actually removed).  Without demodulation the  signal can possibly be spread across multiple DFT bins thus reducing the SNR in the bin that was used to measure the phase. If the signal bandwidth is less than one bin, and the noise is still white, then there would be further increase in SNR by the rate of $\sqrt{N}$ by increasing the number of samples by $N$.
The extent of this spreading due to modulation isn't clear since only the modulation index was provided; we would also need to know the modulation rate and type of waveform in order to predict the occupied bandwidth. With both the modulation index and rate with a continuous modulation waveform the occupied bandwidth can be approximated using Carson's Rule which is applicable to FM and PM given as:
$$BW \approx 2(\Delta F + F_m) \tag{1} \label{1}$$
With:
$BW$: Approximate occupied BW
$\Delta F$: frequency deviation $F_m$: the modulation rate.
And using the modulation index $\beta$ which is the peak phase angle (for sinusoidal PM) and related to $\Delta F$ and $F_m$ according to:
$$\beta = \frac{\Delta F}{F_m} \tag{2} \label{2}$$
From \ref{1} and \ref{2} we can establish the occupied bandwidth from the modulation rate and modulation index directly as follows:
$$BW \approx 2F_m(\beta + 1) \tag{3} \label{3}$$
From this we can confirm if the occupied bandwidth exceeds the resolution bandwidth of the measurement, at which point the SNR from a single bin will be maximized as limited by the DFT noise floor. If the occupied bandwidth is less, and the DFT noise floor still white in close proximity to the signal (meaning the signal is sufficiently stationary over that duration of time), then there will be a further improvement simply by increasing the number of samples.
The resolution bandwidth for each bin in the DFT given in Hz is given by the following formula as the "Equivalent Noise BW":
$$\text{ENBW} = \frac{f_s}{N}w_b \tag{4} \label{4}$$
With:
$\text{ENBW}$: Equivalent Noise Bandwidth (Hz)
$f_s$: Sampling rate (Hz)
$N$: Number of samples in DFT
$w_b$: ENBW of the window function in units of DFT bins
Computing the ENBW of window functions is detailed here: Find the Equivalent Noise Bandwidth.  The equivalent noise BW of the Blackman Harris window is approximately 1.7 bins, so the equivalent noise bandwidth in Hz with the 10 MHz sampling rate and 32768 bins would be:
$$\text{ENBW} = \frac{10e6}{32768}1.7 = 518 \text{ Hz}$$
Using this result in \ref{3} results in the maximum modulation rate of
$$F_m \le \frac{\text{ENBW}}{2(\beta+1)} = 123.3 \text{ Hz}$$
Therefore, if the modulation rate is not significantly higher than 123 Hz, the signal will be entirely within the resolution bandwidth of the DFT and no further phase demodulation would be needed, within the SNR provided by the DFT as was done. If the modulation rate is significantly less, and the signal still stationary over longer duration (eventually phase noise of the signal source will limit this condition), then SNR can be further improved by $\sqrt{N}$ by increasing the number of samples by $N$, reducing the DFT noise floor. If the modulation rate is significantly more than 123 Hz, then there would be additional benefit in SNR to first demodulating the signal to remove the phase modulation. (Also to note that at higher modulation indexes, more signal energy is concentrated in the outer portions of the bandwidth such that there would be a decrease in SNR with further time duration or increasing $N$). When the modulation rate is less than the ENBW of the DFT, even though the signal is phase modulated, given the "courseness" of the DFT, the signal would appears just as a pure tone in this case since the DFT is providing us the average frequency, average phase over the duration of the capture.
So continuing the assumption of $F_m \le 123 \text{ Hz}$, with a phase variations of 0.15° and no motion (meaning this is phase noise); assuming the 0.15° was an rms measurement, that would be equivalent to an SNR as limited by noise due to phase of:
$$-20log_{10}((0.15) 2\pi/360) = 51.6 \text{ dB} $$
This is a 51.6 dB SNR in the resolution bandwidth of 518 Hz; if the signal for which the DFT was computed was real, then the SNR of the signal would be 3 dB higher given it would be at both bins 327 and 32441 and only one of those two bins was used for the measurement. The OP has not stated what bandwidth was used for the original 30 dB SNR metric of the signal, but this approach as outlined can be used to further quantify if the process is noisier than expected, or not. If the variation in phase was anything else such as peak to peak or peak measurements would translate to an even higher SNR, but that should be clarified as well for an accurate assessment. This should also help quantify the SNR requirements for a certain phase accuracy as needed for the final distance measurement.

When such phase noise is small, the small angle approximation $\sin(\theta) \approx \theta$ applies, thus the rms of any phase noise when normalized to the signal (typically given as "dBc") is the SNR when both units are squared to be power terms (hence $20Log_{10}(\theta_{rms}$))


Sampling Dirac function and a DC signal Can we sample the Dirac function?

$$ x(t) = \delta(t) $$

Can we sample a DC signal?

$$ x(t) = 1$$
I think that we can't sample $x(t) = 1$ because the Fourier Transform of $x(t) = 1$ is $2\pi \delta(\omega)$ (the Dirac function)
and we can't sample the Dirac function because it extends to $\infty$ and we cant find $f_\mathrm{max}$.
I find this solution but I am not sure. Can anyone help?

Can we sample the Dirac function?

Strictly speaking: "sampling" would be taking the instantaneous value. Since the Dirac Delta doesn't have a value at $t=0$ (it is not really a function!), NO.
Realistically speaking: an ADC can't measure instantaneous values. That's impossible, because it would require infinite bandwidth of the conversion system (not only of the digital signal), and physical reality doesn't allow for that. So, what a device with limited bandwidth does is measure the (weighted) integral over some (short) time interval. And then, the Dirac isn't any problematic: unlike a value at every point, we know its integral over every interval.
However, when we take physical reality of ADCs into account: there can't be a Dirac Delta in real world, that would require infinite bandwidth... we hence can't measure something that doesn't exist, so yet another: no.

and we can't sample the Dirac function because it extends to ∞ and we cant find fmax.

aaaah! that's not the reason. The signal not being observable is not the same as not being truthfully reconstructable from the digital representation.

Can we sample a DC signal?

yes of course.

I think that we can't sample x(t)=1 because the Fourier Transform of x(t)=1 is 2πδ(ω) (the Dirac function)

No problem with that. Sampling happens in time domain.
It's also not a problem for reconstruction. Nyquist demands limited bandwidth, and your signal very much has limited bandwidth (it has zero bandwidth, to be exact).


Does convolution done via the FFT _undo_ computational gains when signal has zeros every other sample? Suppose we have filter $h[n]$, length $N$. If I use it to do convolution in the time domain, then there will be $N$ multiplications each time during the sliding. (Ignore border conditions for now). However if say $h[n]$ had zeros every other sample, then there would be half that number of multiplications, so ostensibly I save in computational load, since I have half the number of multiplies to do.

So then in conclusion the speed of convolution in the time domain is related to the number of zeros.

Now suppose I did this convolution via the multiplication in the frequency domain using the FFT instead. It seems that I would undo any gains I could have had, zeros or no zeros, because when we do it via the FFT method, it is agnostic to the actual values of the samples of the signal and filter.

Is this a correct assessment? If so, might there then be cases where a good old fashioned convolution in the time domain is in fact much faster than an FFT?

If you know ahead of time that your filter pulse response has zeroes in alternation, then you can use a different FFT routine to speed up the convolution. On the other hand, if you simply pass the filter input and the pulse response to a generic computational subroutine that will spit out the filter output given $x[n]$ and $h[n]$, then yes, the zeroes will be handled no differently than any other inputs and will not save any computation at all.

Example:
Suppose that the (FIR) filter pulse response is $\mathbf h = 1 \quad 0 \quad 1 \quad 0 \quad 1$.  Then, the filter output is

$$y[n] = x[n] + x[n-2] + x[n-4]$$

and so the filter output is alternately
the sum of three successive even-numbered inputs and
the sum of three successive odd-numbered inputs.  So, knowing that
$\mathbf h$ has alternate zeroes, we could proceed as follows.

Create a filter $\hat{\mathbf h}$ by decimating $\mathbf h$
(or downsampling) by a factor of $2$. Thus, $\hat{\mathbf h} = 1 \quad 1 \quad 1$ in the example.
Decimate $\mathbf x$ to create two input streams
$$\mathbf x_e = x[0]\quad x[2]\quad x[4] \quad \cdots ~~\text{and}~~
x_o = x[1]\quad x[3]\quad x[5] \quad \cdots$$
Apply $\hat{\mathbf h}$ to the two streams to get $\mathbf y_e$ and $\mathbf y_o$.
Recombine $\mathbf y_e$ and $\mathbf y_o$ to get $\mathbf y$.

If FFTs are used anywhere in doing the needed convolutions, note that the
FFT of $\hat{\mathbf h}$ needs to be computed only once, and various tricks
such as overlap-and-add etc will have to be applied separately to the two
data streams. Note also that the FFTs can be shorter (though we have to
do two of them), and a crude comparison
says that
$$2\left[\frac{N}{2}\log_2\left(\frac{N}{2}\right)\right] = N\log_2\left(\frac{N}{2}\right) < N\log_2 N.$$
More exact comparisons with actual numbers may reveal a different tradeoff.


What does this paragraph mean? I'm reading this paper about about an algorithm to measure image sharpness, and am confused by these sentences in the 6th paragraph:

It is well known that the attenuation of high-frequency content can lead to an image that appears blurred. One way to measure this effect is to examine the image’s magnitude spectrum $M(f)$, which is known to fall inversely with frequency, i.e. $M(f)\propto f^{-\alpha}$, where $f$ is the frequency and $-\alpha$ is the slope of the line $\log{M}\propto-\alpha \log{f}$.

I don't know what frequency this is referring to as the images this paper and algorithm are concerned with are grayscale, does anyone have a clue?

Detail in images require higher frequency basis functions.  The frequency in this case is measuring fluctuations in intensity as a distance is traversed.  With a lot of detail, a lot of fluctuations, thus higher frequency.
Tamp down the higher frequency and you lose detail, i.e. the image blurs.
The dampening is measured (described) best on a log scale.
Consider this answer Honey Mustard.


Homography Matrix Correlations If the homography matrix of Image 1 to Image 2 is given, and another homography matrix for image 2 to image 3 is given, how do I calculate the homography matrix for image 1 to image 3?
PS : I need the homorgraphy matrices to be 3x3 only because I would be working on a large dataset.

Consider three points $A$ in image 1,$B$ in image 2 and $C$ in image 3.$H1$ is the homography matrix from $A$ to $B$ and $H2$ is the homography matrix from $B$ to $C$.
Then,

$$B = H1A$$

$$C = H2B$$

Substituting for B from first equation

$$C = H2(H1A)$$

Due to matrix associative property,

$$C = (H2H1)A$$

The homography matrix from $A$ to $C$ is given by $H2H1$.


Confusion regarding Nyquist Sampling Theorem The first time Nyquist Theorem was mentioned in class. It stated that we should sample at twice the highest frequency content of the signal. Example: If we wanted to sample $\cos(2 \pi f_0 t)$, the sampling frequency should be at least $2f_0$.

However, in another course. The Nyquist Theorem was stated as such: the sampling frequency should be at least twice the bandwidth of the signal.
Isn't the bandwidth of a single tone cosine 0? which makes the two definition contradictory.

It could be interesting to start from an history of this theorem in Interpolation and Sampling: E. T.  Whittaker, K. Ogura and Their Followers, by Butzer et al.

Putting history aside, the main thing with the basic theorem is that one should refrain to say that the rate  "should" or "has to" be above something. This is could a sufficient, but not a necessary condition. One version is, with an "if": if $X(f)=0$ for $|f| > B$,  sampling at a rate above $2B$ theoretically allows you to recover the signal from the  regularly sampled sequence:

There are milder versions, but this one shows that you ought to be especially careful if you choose, dangerously, $B$ such that $X(B)\ne 0$. For instance with your sine, whose spectrum is not $0$ at $B=f_0$. But is is only an if.

Some non bandlimited signals can still be sampled perfectly, under additional conditions. Some signals can be sampled at $B$ even if they spectrum does not vanish at $B$. And some signals can be sampled at a lower rate. Especially when the signal is bandlimited. @JasonR already pointed to Undersampling, and is dealt with in The theory of bandpass sampling, Vaughan et al., and already discussed in Minimum Sampling Rate of Bandpass signal.

I am not so confident with my understanding of this literature. However, in favorable circumstances with a real bandpass signal, twice the effective bandwidth can suffice, which can be much lower than twice the maximum frequency. And there are other theorems for signals whose spectrum is made of unions of segments of non null frequency.

While mentioning such results, I believe the real issue is different. In some texts, authors assimilate (wrongfully) what they call "the bandwidth" with the distance from the $0$ frequency to the maximum frequency. Especially in texts that do not go into fine details such as the ones stated above and in other answers.

So in your case, I assume that the difference between the two courses can be explained by laziness or mundane talking.

In practice though, you cannot rely solely on this theorem: it requires you to know in advance the spectrum of a continuous signal, which in many cases you have not idea about. Physical modeling  and analog filtering can help you, but in real life, finite precision quantization, signal jitter, noise and especially finite length signals (which cannot have bounded spectra) forster more precaution and a choice of a sampling frequency sufficiently  higher than twice the maximum frequency (or bandwidth if that applies).

A potential additional lecture: Sampling: What Nyquist Didn’t Say, and What to Do About It.


High Level Process, Calculating Cepstrum When You Only Have Spectrum I am not a DSP expert, so I would like someone to check my thought process:

I have a text file representing the spectrum (calculated by some FFT hardware), and I'd like to calculate a cepstrum.  I'm using Octave in the examples below, but it's not required.

My source data is two columns: frequency, amplitude, i.e.:

+----------+----------+
| 0.00E+00 | 6.93E+01 |
| 1.00E+01 | 6.95E+01 |
| 2.00E+01 | 7.38E+01 |
| ...      | ...      |
| 2.00E+03 | 6.51E+01 |
+----------+----------

Load the data in to Octave, and use ifft to put the spectrum into
the time domain.  Then use rceps to calculate the cepstrum:

rceps(ifft(my_data(:,2))); % my_data is the 2000x2 matrix from above
Plot the result of rceps.  The resulting curve should equal cepstrum that would have been produced if I had started with a time domain signal.

Is this a correct way to find the cepstrum when only spectrum (FFT) is available, or, is there a flaw in my thought process?

The cepstrum is defined as:

$C = | \mathcal{F} \log( P(f) ) |^2$

$\mathcal{F}$ is a Fourier transform and $P(f)$ is the power spectrum.  If you've got the power spectrum for the signal there is no need to go back to the time domain.  Just take the logarithm, compute an FFT, take the magnitude squared of the result.

Also, the data you list is strictly positive and real so I'm assuming it's the spectrum.  At the end of your question you put spectrum (FFT).  If what you had was an FFT of a real-valued time domain signal, and not a spectrum, I would expect it to have an imaginary part.


Is this an error in Oppenheim and Schafer's Discrete-Time Signal Processing (3rd edt)? Q: According to Fig 71, y[n]=something+e[n] and hence the transfer function from e[n] to y[n]
denoted as He(z) should be 1. Am I wrong?


Is there any opposite name for Decibel? According to definition of decibel it’s “logarithmic unit of measurement used to express”

So in audio we often convert linear value from range 0 to 1 to logarithmic scale for example like that:

double result = 20.0 * log10(valueInRange_0_1);

(of course we need to menage situation when valueInRange_0_1 = 0, but leave it now, it’s not the case)

So for example for valueInRange_0_1 = 0.5 result will be about -6 and we can call it “-6 decibels”.

But the question is how we can call valueInRange_0_1 = 0.5:

“0.5 OF WHAT” ?

I know in audio we can call it “0.5 gain”, or “0.5 of linear value”. But is there any general term to use as opposite to decibel?

Most of the time, it's the ratio of a physical quantity to a reference. The reference can either be an industry standard (Example: dBSPL, dBm, dBu etc.) or some other suitable quantity (Example: voltage at the input of a filter).

"Physical quantities" in audio are things like pressure, volume velocity, particle velocity, displacement, voltage, current, force, etc. (20 log) or sometimes also power, energy, or intensity (10 log).

For example: a microphone measure "sound pressure" as a physical quantity and to make it logarithmic you use

$$P = 20 \cdot log_{10} (p/p_0)$$

where $p$ is the original sound pressure, $p_0$ a suitable reference pressure, (20 micro Pascals is industry standard) and P the resulting sound pressure level in dBSPL.

Taking a ratio is very important, since it also takes care of the units. In digital signal processing, you often just deal with numbers (without units), but a proper choice of reference is still important. Typically choices are either digital full scale, some suitably calibrated physical quantity or just simply "1" . The reference is important: a dB value without a well defined reference (implicit or explicit) is meaningless.


What image processing operation should I use to remove fuzzy black area I'm performing some rudimentary OCR on some images that I have. The image below is from a VIN tag in a car. The image was originally color with a lot of glare that made some of the characters illegible. I performed the following operations.

Color Invert (tag was black and text was white, I wanted the reverse)
Color Threshold to get rid of the worst of the glare.
Make image black+white only
Erode/Dilate in order to tinker with the text to make it easier for Tesseract OCR to read.
Run Tesseract OCR to extract the characters

The issue that I'm having now is that the black "fuzz" that surrounds the sticker is slowing down Tesseract. What I would like is the name of some operation or a series of steps that I can regularly use in a case like this to get rid of these distortions. I could manually write something to look at squares of the image and if there's way more dark then light, replace with white. There has to be some formal image processing technique to handle this case that I just don't know of.

If I manually remove the black in the margins, my accuracy goes up substantially in Tesseract.

You can perform local averaging/blurring of the picture and then use a threshold to remove all areas that are too dark:
img = cv2.imread('k1aqp.png')
figure(figsize=(20,20))
subplot(2,2,1);
imshow(img)

subplot(2,2,2);
B = cv2.blur(img, (55,55), 0)
imshow(B)

subplot(2,2,3)
mask = B > 150
imshow(mask)

subplot(2,2,4)
img2 = img.copy()
img2[~mask] = 255
imshow(img2)

you need to play around with the blur size and the threshold. But in general, this is one way to go.


Why is the second half of the FFT negative frequencies There are lots of questions here about what the negative frequencies in fft mean, but I'm confused on why the second half of the fft calculation is the correct calculation for a negative frequency.

The equation is
$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2 \pi kn/N}$$

I interpret this as the correlation between a wave with 1 cycle over N, a wave with 2 cycle over N, a wave with 3 cycles over N.

Let say I'm using 10 bins, numpy will tell me the 10th bin is actually representing a negative frequency. But looking at the equation isn't it a correlation between the signal and a signal with 10 cycles over the full period? If they plugged in a negative k into this equation that would make sense. Why does plugging in big values of K give the correlation for a negative frequency?

If you say that plugging in a negative value for $k$ would make sense to you then have a look at this:

$$\begin{align}X[N-1]&=\sum_{n=0}^{N-1}x[n]e^{-j2\pi (N-1)n/N}\\&=\sum_{n=0}^{N-1}x[n]e^{-j2\pi n}e^{j2\pi n/N}\\&=\sum_{n=0}^{N-1}x[n]e^{j2\pi n/N}=X[-1]\end{align}$$

The DFT $X[k]$ is periodic with period $N$, so it doesn't matter if you consider the interval $k\in[0,N-1]$, or $k\in[-N/2+1,N/2]$ (assuming $N$ is even).


How to derive the $\log$ of fraction of two Gaussian distribution I am designing a function $A(x)$ that map $x\to \{0,1\}$. The function $A$ can be expressed as
$$A(x)=\begin{cases}
1 \text{  if $p_1(x) \ge p_2(x)$}\\
0 \text{  otherwises}
\end{cases}
$$
where $p_1$ and $p_2$ are Gaussian distribution:

$p_1(x)=\frac{1}{\sqrt{2\pi}\sigma_1}\exp\left( {-\frac{(x-\mu_1)^2}{2\sigma_1^2}} \right)$;

$p_2(x)=\frac{1}{\sqrt{2\pi}\sigma_2}\exp\left( {-\frac{(x-\mu_2)^2}{2\sigma_2^2}} \right)$.

Could you help me to derive a short form of $A(x)$ ?

This is my solution. However, it maybe wrong

I derived the bellow eq.
$$A(x)=\begin{cases}
1 \text{  if $\log \frac{p_1(x)}{p_2(x)}>0$}\\
0 \text{  otherwises}
\end{cases}
$$
where $\log(e)=1$

We have
$$\log \frac{p_1(x)}{p_2(x)}=\log \frac{\sigma_2}{\sigma_1}+\left( {\frac{x-\mu_2}{\sqrt{2}\sigma_2}} \right)^2-\left( {\frac{x-\mu_1}{\sqrt{2}\sigma_1}} \right)^2$$

Hence, my final solution is
$$A(x)=\begin{cases}
1 \text{  if $\left (\log \frac{\sigma_2}{\sigma_1}+\left( {\frac{x-\mu_2}{\sqrt{2}\sigma_2}} \right)^2-\left( {\frac{x-\mu_1}{\sqrt{2}\sigma_1}} \right)^2 \right) \ge 0$}\\
0 \text{  otherwises}
\end{cases}
$$

Thanks in advance

Your final condition is correct, but it doesn't show explicitly for which values of $x$ it is satisfied. This can be worked out as follows. From what you have it follows that $p_1(x)\ge p_2(x)$ is equivalent to

$$0\ge x^2(\sigma_2^2-\sigma_1^2)-2x(\mu_1\sigma_2^2-\mu_2\sigma_1^2)+\mu_1^2\sigma_2^2-\mu_2^2\sigma_1^2-2\sigma_1^2\sigma_2^2\log\left(\frac{\sigma_2}{\sigma_1}\right)\tag{1}$$

Let's first consider the case $\sigma_1=\sigma_2$. From $(1)$ it follows that in this case the condition on $x$ becomes

$$x\ge \frac{\mu_1^2-\mu_2^2}{2(\mu_1-\mu_2)}
= \frac{\mu_1+\mu_2}{2}\tag{2}$$

where we assume $\mu_1\neq\mu_2$ because otherwise $p_1(x)=p_2(x)$, which makes no sense.

In the case $\sigma_1\neq\sigma_2$, we have to consider $(1)$ as a quadratic inequality in $x$, which is satisfied with equality for the following values of $x$:

$$\begin{align}x_1&=\frac{a}{2}-\sqrt{\frac{a^2}{4}-b}\\
x_2&=\frac{a}{2}+\sqrt{\frac{a^2}{4}-b}\end{align}\tag{3}$$

with

$$\begin{align}a&=\frac{2(\mu_1\sigma_2^2-\mu_2\sigma_1^2)}{\sigma_2^2-\sigma_1^2}\\
b&=\frac{\mu_1^2\sigma_2^2-\mu_2^2\sigma_1^2-2\sigma_1^2\sigma_2^2\log\left(\frac{\sigma_2}{\sigma_1}\right)}{\sigma_2^2-\sigma_1^2}\end{align}\tag{4}$$

Depending on the sign of $\sigma_2-\sigma_1$ we have two conditions equivalent to $p_1(x)\ge p_2(x)$:

$$\begin{align}\sigma_2>\sigma_1:\quad & x_1\le x\le x_2\\
\sigma_2<\sigma_1:\quad & x\le x_1\text{ or }x\ge x_2\end{align}\tag{5}$$

The inequalities in $(5)$ state explicit conditions on the variable $x$ that can be easily evaluated.

EDIT: As suggested by Dilip Sarwate in a comment it can be instructive to distinguish another special case: $\sigma_1\neq\sigma_2$ and $\mu_1=\mu_2=\mu$. In that case the limits $x_1$ and $x_2$ of the decision regions given by $(3)$ lie symmetrical to the common mean $\mu$, and $(3)$ simplifies to

$$x_{1,2}=\mu\pm\sqrt{\frac{2\sigma_1^2\sigma_2^2\log\left(\frac{\sigma_2}{\sigma_1}\right)}{\sigma_2^2-\sigma_1^2}}$$

Of course, the conditions $(5)$ remain valid.


Speech Compression - In LPC how does the linear predictive filter work on a general level? Hi I'm taking a multimedia systems course and I'm preparing for my exam on tuesday.  I'm trying to get my head around LPC compression on a general level, but I'm having trouble with what is going on with the linear predictive filter part.  This is my understanding so far:
LPC works by digitising the analog signal and splitting it into segments.  For each of the segments we determine the key features of the signal and try to encode these as accurately as possible.  The key features are the pitch of the signal (i.e. the basic formant frequency), the loudness of the signal and whether the sound is voiced or unvoiced.  Parameters called vocal tract excitation parameters are also determined which are used in the vocal tract model to better model the state of the vocal tract which generated the sound.  This data is passed over the network and decoded at the receiver.  The pitch of the signal is used as input to either a voiced or unvoiced synthesiser, and the loudness data is used to boost the amplitude of this resulting signal.  Finally the vocal cord model filters this sound by applying the LPC coefficients which were sent over the network.
In my notes it says that the vocal tract model uses a linear predictive filter and that the nth sample is a linear combination of the previous p samples plus an error term, which comes from the synthesiser.

does this mean that we keep a running average of last p samples at both the encoder and decoder?  So that at the encoder we only transmit data that corresponds to the difference between this average and actual signal?

Why is it an a linear combination of these previous samples?  My understanding is that we extract the loudness, frequency and voiced/unvoiced nature of the sound and then generate these vocal tract excitation parameters by choosing them so that the difference between the actual signal and the predicted signal is as small as possible.  Surely an AVERAGE of these previous samples would be a better indication of the next sample?

If there are any holes in my understanding if you could point them out that would be great!  Thanks in advance!

LPC voice coders (starting with the old LPC10 standard, which seems to be the one you refer to here) are based on the source-filter model of speech production. Speech can be characterized by the following properties:

The raw sound emitted by the larynx (through vibration of the vocal folds, or just air flowing through it, the vocal folds being opened).
The transfer function of the filter achieved by the articulatory system, further filtering this raw sound.

Early LPC coders (LPC10) adopt the following model of those two steps:

The larynx emits either white noise, characterized by an amplitude $\sigma$ ; or a periodic impulsion train, characterized by an amplitude $\sigma$ and a frequency $f_0$
The transfer function of the articulatory system is of the form $\frac{1}{1 - \sum_k a_k z^{-k}}$, and is thus entirely characterized by the coefficients $a_k$

The principle of early LPC coders (such as LPC10) is thus to estimate these parameters from chunks of the incoming audio signal ; transmit these over the network ; and have a sound generator reproduce a sound from these parameters on the receiver. Observe that in this process, nothing of the original audio samples is actually transmitted. To make a musical analogy, it's like listening to a piano performance, transcribing it, sending over the score, and getting someone to play it on the other end... The result on the other end will be close to the original performance, but only a representation has been sent over.

  does this mean that we keep a running average of last p samples at
  both the encoder and decoder?

No, this is not how it works. On the encoder side, we run a process known as AR (autoregressive) estimation, to estimate the set of coefficients of the AR filter that matches best the spectral envelope of the input signal. This is an attempt at recovering the coefficients of the filtering that has been performed by the articulatory system.

These coefficients are sent over the network (along with pitch, the voice/unvoiced flag, and loudness). The decoder uses those coefficients to filter a synthetic excitation signal, which can be either white noise (unvoiced frame) or a comb of periodic impulsions (voiced frame). These coefficients are used in the following way on the decoder to recover the output signal $y(n)$:

$$y(n) = excitation(n) - \sum_k a_k y(n - k)$$

Observe that since this is an all-pole, IIR filter, the samples linearly combined are previous samples produced by the filter.

  So that at the encoder we only transmit
  data that corresponds to the difference between this average and
  actual signal?

There is no "averaging" and there is no transmission of a difference signal. The decoder is a sound synthesizer with some parameters ; and the encoder searches for the set of parameters for this synthesizer that most closely matches the input signal.

  Why is it an a linear combination of these previous samples?

Other options would be possible indeed here, but the advantage of using an auto-regressive model are the following:

In this model we use an all-pole IIR filter to simulate the articulatory system. This is a good model, since it can capture strong spectral dips and peaks with a small number of coefficients ; and the articulatory system is indeed capable of emphasizing / attenuating narrow bands of frequencies (See formants). If we had used an all-zeros FIR filter, we would have needed many more coefficients to accurately capture the kind of filter responses that the articulatory system can achieve.
On the encoder end, the problem of of estimating the filter coefficients is computationally efficient - it can be done using the Levinson-Durbin recursion on the first few samples of the autocorrelation. More complex linear models (ARMA) or non-linear models are more computationally expensive or intractable.
On the decoder end, the synthesis is very simple indeed - we just need enough memory to keep track of the n previous samples emitted by the decoder.

  Surely an AVERAGE of these previous samples would be a better indication of the next sample?

This is not true. For example, assume the input signal is a sine wave. The prediction $\hat{y}(n) = (2-\omega^2) y(n - 1) - y(n - 2)$ has zero error, while $\hat{y}(n) = \frac{y(n - 1) + y(n - 2)}{2}$ is wrong most of the time (in particular, it does not capture the fact that the sine wave is strictly increasing half of the time, strictly decreasing the rest of the time, so $y(n)$ shouldn't be "between" the previous values). Speech signals are not sinewaves of course - but you can squeeze much more modelling power out of a model with $p$ parameters (order p AR), than out of a model with zero (averaging).

Also, it's worth keeping in mind that there's a mathematical definition of "better" (minimizing the expected value of the squared of the error), which yields a procedure for finding the optimal value of the coefficients.