How to increase the spectral resolution? My question is about the spectral resolution of a discrete signal. Each sample of my signal is made up with 2^n frames sampled at 44.1 kHz.

So, when I want to know the spectral resolution, I calculate : 44100/number_of_frames. With 2048 frames, my spectral resolution is around 20Hz. But, when I take a look to the bands filtered by an equalizer, the band with in the low frequencies are around only 5 Hz (16 Hz -> 20 Hz -> 25 Hz -> 32 Hz...). How is it possible ?

I thought about using "zero padding", but, even if it will help me to have a better location of each peak of the the spectrum analyze, this method don't magically increase the spectral resolution.

I also thought about increasing the number of frames analyzed. But, to get a spectral resolution of 5 Hz with a signal sampled at 44.1 kHz, I would need 8192 frames and it would represent 185 ms. It's very far from a pseudo real-time analysis and a singer who would listen his voice after this analyze while he is singing would hear this "delay".

So, what is the solution ?

Thank you for all your reply

You seem to have a good grasp of the tradeoffs. When using short-time Fourier analysis like you are, there is a version of the uncertainty principle at play. Increasing your time resolution (in your case, using a shorter DFT) results in coarser frequency resolution, and vice versa. That is, the time-bandwidth product is a constant.

The way to increase your STFT's spectral resolution is to increase the duration of time that the transform covers, as you noted. If you truly need to be able to resolve frequencies that precisely (within a few Hz of one another), then you need to observe them for a long enough period of time to discern them. If you know a priori some characteristics of your signal, and conditions are favorable (i.e. SNR is high enough), then you might be able to get the job done with coarser resolution (and therefore a shorter transform).

For instance, if you know that your signal is likely to be a single tone somewhere in a particular band, and you want to know its frequency precisely, then you don't necessarily need to use a really long DFT. Instead, you can use a shorter DFT, then use peak interpolation techniques to give a sub-bin estimate of where the peak actually lies.


Estimating onset time of a tone burst in noise? What techniques might one use to estimate the onset time of a sinusoidal tone burst in a noisy signal?

Assume the tone burst has a known fixed frequency (but unknown phase) and a very sharp rise time, and that the goal is to estimate the onset time within better than half the rise time, and/or one period of the frequency of the tone, if possible.  How might the estimation techniques change if the S/N ratio is very low (much less than 1)?

Added: Assume the tone burst is of unknown length, but longer than a small multiple of the rise time and the frequency period.

Added: A DFT/FFT shows the very probable existence of a tone.  The problem is figuring out exactly precisely where in the FFT window the tone (or perhaps multiple tone bursts of the same frequency) may have started within the FFT window, or determining if the current tone started outside that DFT window, provided I have all that additional time domain data.

Radar pulse detection accuracy is closer to the resolution I need, except I only have an edge, as the tone is of unknown length, and, other than a known rise time, unmodulated.  Narrow band pass filters distort the rise time, and thus kill edge arrival estimation resolution.

As we've been discussing in the comments, the Goertzel algorithm is the usual way to detect a tone in noise.  After the discussion, I'm not sure it's quite what you are after (you want the onset time), but there seemed to be confusion over how the Goertzel algorithm might be applied to your problem, so I thought I'd write it up here.

Goertzel Algorithm

The Goertzel algorithm is good to use if you know the frequency of the tone you are looking for (call it $f_g$), and if you have a reasonable idea of the noise level so that you can select an appropriate detection threshold.

The Goertzel algorithm can be thought of as always calculating the output of ONE FFT bin:

$$y(n) = e^{\jmath 2\pi f_g n} \sum_{k=0}^n x(n) e^{-\jmath 2\pi f_g k}$$

where $f_g$ is the frequency you're looking for.

The Wikipedia page has a better way to calculate this.

Here is a (feeble) Scilab attempt at implementing it:

function [y,resultr,resulti] = goertzel(f_goertzel,x)
realW = 2.0*cos(2.0*%pi*f_goertzel);
imagW = sin(2.0*%pi*f_goertzel);

d1 = 0;
d2 = 0;

for n = 0:length(x)-1,
    y(n+1) = x(n+1) + realW*d1 - d2;
    d2 = d1;
    d1 = y(n+1);
    resultr(n+1) = 0.5*realW*d1 - d2;
    resulti(n+1) = imagW*d1;
end
endfunction

Consider the signal with $f = 0.0239074$ and $\phi = 4.4318752$ :

$$x = \sin(2\pi fn + \phi) + \epsilon(n)$$

where $\epsilon(n)$ is zero-mean, unit variance Gaussian white noise.

In this example, the tone starts one third of the way into the signal at index 1001.

If we run the Goertzel algorithm on it with $f_g = f - 0.001$ then we get the top two traces of the figure.

If we run the Goertzel algorithm on it with $f_g = f$ then we get the bottom two traces of the figure.

The four traces are:

$x$ (blue) and $y$ (red) for $f_g = 0.0229074$
The resulting $\sqrt{resultr^2 + resulti^2}$
$x$ (blue) and $y$ (red) for $f_g = 0.0239074$
The resulting $\sqrt{resultr^2 + resulti^2}$ (solid line) and the first result (dashed line).

As you can see, the case where the tone we are interested in is present peaks at about 250.  If we set the detection threshold at about half this value (125), then the detection occurs (the square-rooted value is greater than 125) at about index 1450 --- 450 samples after the tone started.

This threshold (125) will not cause a detection in the other case (for this run, anyway), but the maximum value of that output is 115.24, we we cannot reduce the threshold too much without getting a false detection.

Reducing the threshold to 116 will cause detection in the true case (for this run) at index 1401... but we run the risk of more false alarms.


How to do Frequency Scaling on an Audio File Please excuse me if my terminology is wrong. I'm from a music production background and have no experience in signal processing.
I was wondering if it was possible to stretch out the overtones (partial tones) of sound through digital signal processing.
In the natural harmonic series  or the overtone series, there is the fundamental frequency. Then you have overtones at integer multiples of the fundamental frequency.
Is there a software which would let me alter the ratio between the different overtones of the sound?
For example, If I have an audio file of 0-20,000 hertz. How can I take the first 10,000 hertz of the audio spectrum and stretch it out to 20,000 hertz so that all the overtones are spaced wider and the ratio is doubled.
It's called Frequency Scaling or Pitch Scaling I believe. I do understand the difference between Frequency shifting and Frequency Shifting and I'm curious to experiment with Frequency Scaling and not Frequency shifting.
I'm looking to do something like this.
Thank you.

Everything is possible with advanced technology, due diligence and perseverance. To discuss the means to reach your goals requires even more detailed description of target scenarios.
If you have a pure tone -- I mean, a musically pure single note performed by a good musician with a quality instrument or an opera singer -- you can analyze this sound by spectral methods and even real-time-process this sound in order to shift its frequency or change its overtones any way you like. Then, you can implement these transformations with a digital signal processor or conventional computer.
The task is more difficult, if you have a musical piece recorded and want to process individual parts of the performance through submitting them to a "pitch scaling" as defined in your Correlogram Museum reference. Even to "pitch scale" a recorded or live-performed triad, not a single note, is an algorithmic and methodological challenge.
In the Correlogram Museum examples, the performer plays on a synthesized piano-like instrument. The DSP unit most certainly takes part in producing the sounds of this instrument, but the pipeline of sound production is inverse w.r.t. your scenario: the sound is generated with "scaled" overtone frequencies, whereas you have a sound with overtones prepared early and want to "post-process" the sound in order to change its overtone content.
I admit there are design projects pertinent to the scenario of your question. Ask colleagues, browse the literature, search the Internet.


Practical insightful time domain functions This question might be a little bit different, but i am having hard time to find some sort of a list/book/website/advice, with very practical time domain functions that provides useful insgihts for a series of numbers in time domain (frequencies of less than 1Hz)

Examples - FFT to find periods in a noisy signal, Autocorrelation  for finding an interesting correlations between 2 curves.

Very curious to know about more, maybe more functionality about prediction and information that is useful and critical in real life.

I've written a series of three blog articles about three families of time domain frequency calculation for single pure tones, real or complex.  You can find them here:

Exact Near Instantaneous Frequency Formulas Best at Peaks (Part 1)
Exact Near Instantaneous Frequency Formulas Best at Peaks (Part 2)
Exact Near Instantaneous Frequency Formulas Best at Zero Crossings

The applicability in DSP is rather limited, but if this is what you are doing, these are the formulas you want.


Can we use ICA (Independent component analysis) on a noisy signal If a receiver received a noisy signal can we use ICA to separate the noise from the signal ?

Yes and no.

The original ICA methods were not made for the case where there is noise on the signal. However, there exist relatively straightforward extensions that work for noisy signals (e.g. [*]). In this case, you can expect some degree of noise suppression but it won't be perfect.

That means that ideally, your desired signal components will stick out enough from the noise for the ICA to recover them. Then again, you cannot expect this recovery to be perfect, there will be some estimation errors. The more noise, the worse the quality of your signal estimate. If you have weak signal components you should be prepared for the chance that they will become invisible under the noise floor. Give it a try I would say.

[*] Hyvarinen, Aapo. "Gaussian moments for noisy independent component analysis." IEEE signal processing letters 6.6 (1999): 145-147.


Convolution between two vectors. Length and normalization I have an RIR vector $h[n]$ with $N$ samples and an audio source $x[n]$ with $M$ samples. I wish to simulate a 5 seconds audio segment with $x[n]$ randomly located within (timewise).

Using MATLABs conv(x,h) I am getting a result with vales in the range $[-0.3852,0.3242]$.

Using pythons np.convolve(x,h) I get a result with vales in the range $[-12621.9,10624.08]$, which also sounds bad on the headset (I am assuming due to cutoffs).

I do not know where is the difference comming from, as both $h$ and $x$ are the same before the convolution. Normalizing the output of the python version by:

output=np.convolve(x,h)
output=output/len(output)

fixes the values. This is true for normalizing by either len(output), len(x) or len(h).

Now I am confused about the best method of action.

For a 5-second segment recording, do I have to generate a 5 - second length $h$?
Is it best to first pad both $h$ and $x$ with zeros and then convolve or should I convolve and than allocate randomly within a 5 seconds zeros vector?
Is it at all reasonable to normalize here?
with respect to the former 3 questions, do I normalize by len(output), len(x), len(h) or the number of samples within a 5 seconds segment?

I am aware that there may be more than one correct answer. I am looking for the pros and cons of each course of action and what is the best way to achieve my target.

The scaling difference between your two output is exactly $32767= 2^{15}-1$ which is exactly the maximum amplitude of a 16-bit signed integer. I'm guessing, the difference is how you import the audio into the program: Matlab's $audioread()$ typically normalize the data to $[-1,1]$, i.e. divide by 32767. It would seem that whatever Python method you are using doesn't do this.

Normalization has nothing to do with the length of the filters and it isn't affected by any zero padding. It's all about the scaling conventions of your inputs and outputs, how do you interface with drivers and/or files.


Why do wavelets integrate to 0 and how do they match a signal? I have been reading about the Wavelet transform recently and its relationship to the Fourier transform.  From what I understand the wavelet transform represents signal data with many short-lived functions (wavelets).  The major difference between the 2 is that the Fourier does not sample in time, but only in frequency.
 I have  couple questions about the wavelet method:

1) Why is it required that wavelet functions integrate to 0?

From my understanding, wavelets are shifted along the time domain and compressed or stretched to match the frequency of the signal at that point in time.  So wavelets are able to capture frequency at different places in time.

2) How does the wavelet capture peaks that are higher or lower in the original domain?  For example see the image below that I have simulated in R.  The dashed line represents the wavelet transformed approximation.  How  is this dashed line able to capture peaks that are higher or lower if it is only shifted and compressed or stretched in the x direction?

For the first question, in words. In the continuous wavelet context, you want to represent a class of signals as a weighted sum of translations and scaling of a single wavelet shape, without losing information. This means that these coefficients can be turned back into the original signal exactly, by a combination of "another" wavelet. For different reasons, if we want this other wavelet to be the initial one, this self-invertibility imposes some restrictions on the possible wavelet candidates. One is the so-called "admissibility" condition, that imposes a zero-mean wavelet (with as well a fast decay for high frequencies). When the wavelet is discretized, this further reduces the diversity of potential wavelet shapes, and still inherit from the zero-mean, or integration to zero. You can look at my other answer here: Admissibility condition for wavelets.

In practice, for discrete or limited-length data, one can tolerate that linear wavelet decompositions are made when the  integral is small enough. Moreover, some multiscale schemes loosely called wavelets (some are non linear) don't even bother. The low-pass, non-zero mean part of data is generally cast to a so-called scaling functions or signal approximation.

For the second question, shift and scale are inherent to the wavelet frame structure. If the signal is multiplied by a constant, a linear wavelet decomposition will multiply its coefficients by the very same constant, as you correctly guessed.


Removing baseline drift from ECG signal I am trying to design a high pass filter to remove baseline drift from an ECG signal.
the baseline drift is of very low frequency like 0.3Hz or so with an amplitude of 25% of the ECG signal. But the bandwidth of ECG signal itself is 0.5Hz to 150Hz. I have been trying fdatoolbox in matlab to design the HPF but I was not able to remove the baseline drift accurately.Please help on how can i achieve a low order filter to eliminate baseline drift.
I have got it using higher order filters of order above 1000. But I'm looking for a way if I can get lower order filter,like order of less than 30, to do the same.Thanks for your help..

You usually use a DC notch filter. This is a recursive filter that should kill a very narrow band of frequencies around DC. One possible implementation is a first order DC notch with transfer function

$$\frac{1 - z^{-1}}{1 - \lambda z^{-1}},$$

where $\lambda$ is some number very close to 1, say 0.99. A to demonstrate this in MATLAB,

b = [1 -1];
a = [1 -0.99];
freqz(b,a);

Notice that the phase response around DC gets pretty non-linear. Such is the nature of these recursive filters.

The closer you make $\lambda$ to 1, the narrower the notch gets and the more distorted your phase response becomes around DC.

For ECG signals, this will usually do. Experiment with $\lambda$ values that will suit your needs. Sometimes, if you make the notch too narrow, you won't get all off the drift, so 0.99 is usually a good starting point.


DC component of a discrete filter I know that many books and papers talk about the DC offset/DC component of a filter. How do we define the DC offset mathematically, for the case of discrete filters?

Here is a (working) link to a paper relevant to the OP: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.80.2334&rep=rep1&type=pdf

It looks like when they say zero DC, they mean "band pass filter" (see page 4). And when they say DC component they just mean frequency response at zero frequency. For instance, see Eq. (20) where they evaluate $G_b(0)$. That's the answer to your question: you can define the DC offset of the filter response as its gain at zero frequency.


EQ - Calculate the power per Octave of a signal through FFT I'm working on a audio Equalizer, and I'm not sure if the way I obtain the power per octave of my signal is good.

Here is the code I have :
(I use the third octave bands, which is initialized earlier in the program)

void    Record::ProcessData(short* Buffer, size_t BufferSize)
{
    /* BufferSize is usualy equal to 131072 */
    /* Initialize the signal and the Hanning Window */
    Aquila::SignalSource signal = Aquila::SignalSource(Buffer, BufferSize, 44100);
    Aquila::HannWindow HannWindow(BufferSize);

    /* FFT */
    /* The library overload the + operator, so the signal is multiplied with the window */
    auto fft = Aquila::FftFactory::getFft(BufferSize);
    Aquila::SpectrumType spectrum = fft->fft((signal + HannWindow).toArray());

    /* Set the imaginary part to 0 */
    for (size_t i = 0; i < BufferSize; i++)
        spectrum[i]._Val[1] = 0;

    /* OctaveRange is a class with the minimum, center and maximum frequency of an octave band */
    list<OctaveRange>::iterator it;
    map<double, double> sumMap;
    double db, f, key;

    /* Initialize the map values to 0 */
    for (it = thirdOctave.begin(); it != thirdOctave.end(); it++)
        sumMap[(*it).getCtr()] = 0.0;

    /* The calcul of the power per octave starts from here */
    for (size_t i = 0; i < (BufferSize / 2); i++)
    {
        /* f = frequency at bin i */
        f = (double(i * 44100) / BufferSize);
        for (it = thirdOctave.begin(); it != thirdOctave.end(); it++)
        {
            key = (*it).getCtr();
            if (f >= (*it).getMin() && f < (*it).getMax())
                sumMap[key] += pow(spectrum[i].real(), 2) + pow(spectrum[i].imag(), 2);
        }
    }

    for (it = thirdOctave.begin(); it != thirdOctave.end(); it++)
    {
        key = (*it).getCtr();
        /* db is the power (in dB) per third octave. */
        db = 20.0 * std::log10(sqrt(sumMap[key]);
    }
}

When I record pink noise with this I obtain some believable values, but I'm still not sure if this is because the code is good, or if it's a happy coincidence.

So my question is simple :
Is this the good way to obtain the power per octave band of a signal ?

Feel free to ask if something is not clear enough.
Thanks for your time.

Difficult to say without a reference to the library functions and classes that you are using. In general it looks okay (if fairly inefficient).

A few things to look at:

The Hamming window needs to be multiplied with the signal, not added to it. Whether the code is correct really depends on how the "+" operator for the clasess is defined
db() can be 10*log or 20*log, whether you need the sqrt depends again on how its's defined in the library
Your octave band filters are bascially "brick wall" filters. Most spectrum analyzers do have a finite slope filters so any one frequency contributes to multiple bands (just by different amounts). Whether that's a problem or not depends on the application

Re point 3: consider the output of the analyzer when you input a sine wave and slowly sweep the frequency from 1kHz to 2 kHz. With a finite slope filter you will see the energy in the 1 kHz octave band slowly going down and the one in the 2 kHz band slowly coming up. They will be equal when the frequency is 1414 Hz. In your case the 1 kHz octave band energy will be constant and once you reach 1414 Hz all the energy will immediately flip to the 2 kHz band.


Contours matching - finding contours displacement I found contours on two images with same object and I want to find displacement and rotation of this object. I've tried with rotated bounding boxes of this contours and then its angles and center points but rotations of bounding boxes don't tell about contour rotation correctly because it's the same for angles a+0, a+90, a+180 etc. degrees. Is it any other good way to find rotation and displacement of contours? Maybe some use of convex hull, convexity defects? I've read in Learning OpenCv about matching contours but it hasn't helped. Could someone give some example?

Examples:

I want to detect for example pink square, and in the second case pen. Other examples could be squares with some holes, stars etc. As I said I want to make some uniwersal thing. Any suggestions are appreciated because I want to test as many methods as it's possible.

Do you have to worry about a difference in scale between contours?  If not than you can simply find the centroid of each contour, and compute displacement by subtracting one from the other.  Then you can calculate the principal axes of the contours, and find the rotation angle between them.

If scaling is involved, then you can calculate the scale factor by taking the ratio of the corresponding principal axes.


Do we use closed loop or open loop information in Bode plot, Nyquist plot and Root Locus? When analyzing a control system via Bode, Nyquist and Root-Locus, are we using open or closed loop information? I'm not able to understand when the poles of the open loop are used and when the one of the closed loop are. I think that what interests us is the information about closed loop, as that determines whether the final system is stable or not, but I'm pretty confused with this issue.

The bode plot is just a plot showing the frequency representation of your system. Any transfer function has one, hence you can use it to see the response of both closed or open loop. Nyquist and the root locus are mainly used to see the properties of the closed loop system. The root locus shows the position of the poles of the c.l. system as the gain of your controller changes. While nyquist diagram contains the same information of the bode plot. But instead of splotting two real graphs of the amplitude and the phase of the complex function the nyquist plot shows the same function in the complex plane


Range and motion in received signals by radar I've got two signals which are consecutive received signals by a radar system. The signals' values are complex (in form $x+\textrm iy$).
I have generated a trasmitted signal using the info available.

Now,  I have to calculate the range and motion(if any) of the objects using comparision of these signals in MATLAB. Comparing the real part using xcorr in MATLAB, I got the range of the objects.
I am strugling with motion of the objects. How can I find if the objects found (4 in my case) were moving or not. I think it has to do something with the phase comparision of the signals.

Any help is appreciated.

You will need more than 2 received pulses to get the Doppler between your sources. What is your ADC sampling rate and PRF?


Correct transfer function with 1st order IIR filter What diffreneces it makes in magnitude and/or phase response when using

$$H(z) = \frac{b_0 + b_1z^{-1} + b_2z^{-2}}{a_0 + a_1z^{-1} + a_2z^{-2}}$$ by adding $b_2=0$, $a_2=0$

instead of

$$H(z) = \frac{b_0 + b_1z^{-1}}{a_0 + a_1z^{-1}}$$

for 1st order filter?

All software I'm using for to build VST plug-ins (SynthEdit, FlowStone, etc.) implements a biquad (BLT) module for coefficient input and now it looks like I can't get 90° phase response with a filter designed using Octave (by Octave plot phase is 90°).

EDIT: Ok, digged the web a bit more and found these three (1 2 3) articles published at analog.com. Is the situation same with one pole allpass filter vs biquad allpass filter (what are the phase responses for these?)?

If the formula becomes the same, the response becomes the same – there's no difference to be expected.


Apply the Discrete Fourier Transform as a Matrix Multiplication in MATLAB 0
I have a vector x of length N x 1, I need to perform the iDCT operation for it using MATALB. I performed that using the below method:
x = [1+i;1-i;-1+i;1+i];    %Assume N = 4
y = sqrt(4)*idct(x);

Then I performed it using the other way such that:
x = [1+i;1-i;-1+i;1+i];    %Assume N = 4
D= dctmtx(4)/sqrt(4);
y = D'*x;

I don't find the value of y equals in the two ways ?? However for the case of of FFT operations, that gives the same results!
Is there a special case for DCT operation?

Pay attention that by default MATLAB use DCT Type II hence the inverse is basically DCT Type III:
vX = [1 + 1i; 1 - 1i; -1 + 1i; 1 + 1i];    %Assume N = 4
vY = dct(vX);
mD = dctmtx(length(vX));
vYY = mD * vX;

vYY ./ vY
max(abs(vY - vYY))

vY = idct(vX);
vYY = mD.' * vX;

vYY ./ vY
max(abs(vY - vYY))

The result:
ans =

   1.0000 + 0.0000i
   1.0000 + 0.0000i
   1.0000 - 0.0000i
   1.0000 + 0.0000i

ans =

   4.9651e-16

ans =

   1.0000 - 0.0000i
   1.0000 - 0.0000i
   1.0000 + 0.0000i
   1.0000 + 0.0000i

ans =

   4.5776e-16

So everything works as expected.
Pay attention that in your code using the transform you multiply by factor sqrt(4) and for the matrix form you divide by this factor.
So:
vY = sqrt(4) * idct(vX);
vYY = (mD' / sqrt(4)) * vX;

vY ./ vYY

ans =

   4.0000 + 0.0000i
   4.0000 + 0.0000i
   4.0000 + 0.0000i
   4.0000 - 0.0000i

So just don't divide and use the same factor and all is good.


Convolution Based Image Interpolation From a book, I learned that image interpolation (or upsampling) can be written as a convolution like this:
$$g(i,j) = \sum_{(x,y)}f(x,y)h(i-rx,j-ry)$$
while $r$ is the upsampling rate.

But I have a time understanding this formula, since AFAIK, the basic idea of image interpolation is to bring near-by pixels into consideration, so I thought a convolution like below might suffice:
$$g(i,j) = \sum_{(x,y)}f(x,y)h(i-x,j-y)$$
What is the point of $r$ in the formula?

Image interpolation is not (exactly) the same as upsampling. Suppose that you have a grid $\Lambda$ which is a $640 \times 480$ grid of pixels. Each pixel has a location $(x,y) \in \Lambda$. Each pixel also has an intensity $I(x,y)$. This could be a gray-scale value or an RGB triplet. Now suppose that you want to know the intensity of a location that is not centered on one of the pixels. For example, suppose that you wanted the intensity at location $(17.35, 188.42)$. In this case, you would use the neighboring integer pixel intensities in an appropriate interpolation scheme (nearest neighbor, bilinear, bicubic, Lanczos etc.) to get the sub-pixel intensity. This is image interpolation.

Image upsampling is when you want to go from a grid $\Lambda$ that is $640 \times 480$ to a larger grid $\Gamma$ that is for example, $1280 \times 960$ pixels. In such a case, a simple algorithm is create the new grid filled with zeros. Then, if there is a pixel value in the original image at $(x,y)$, then take its intensity I(x,y) and fill it in the upsampled image at $I'(2x,2y)$. The $2$ is the $r$ (rate) referred in your formula. Then, we use an appropriate kernel (the $h$ in your formula) to filter the image and complete the upsampling. NOTE: This is a simple explanation. Please see Gozalez, Woods and Eddins, Digital Image Processing for a better treatment.


Reducing sampling rate by a non-integer factor I have data that is sampled at 12 kHz. The downsample function in software such as MATLAB only allow you to downsample by an integer factor, i.e. from $12\to{6,4,3,2,1,\ldots}$ kHz. How would I downsample by a non-integer factor, say, from $12\to{10,9,8,7,5,...}$ kHz (pseudocode/description is fine)? You can assume that the signal has been appropriately lowpass filtered prior to downsampling. Also, the desired new sampling rate (after downsampling) is a "nice" number (i.e., I'm not trying to do $12 \to 9.7432$ or something like that).

I can recommend the resample command in MATLAB. It's a combination of up- and down-sampling. To get from 12kHz to a sampling rate of 9kHz, you upsample by 3 and downsample by 4. resample just needs your signal and the two factors for up- and down-sampling.  Remember, when downsampling, to first low-pass filter out everything between 4.5 kHz and 6 kHz. Very little energy should exist above the resulting Nyquist frequency.


Get the discrete-time poles and zeros from continuous-time poles and zeros How do you implement the following function:

$$[Z^d, P^d, K^d] = \text{fcn} \,(Z^c, P^c, K^c),$$

where $Z^c = [z^c_m, ..., z^c_1]$, $P^c = [p^c_m, ..., p^c_1]$, and $K^c$ are zeros, poles, and gain of the continuous-time TF, and similarly $Z^d$, $P^d$, and $K^d$ are zeros, poles, and gain of the equivalent discrete-time TF with sampling time of $T_s$.

I have tried the following but it did not give me the correct answer:

$Z^d = e^{Z^c \times T_s}$

$P^d = e^{P^c \times T_s}$

$K^d = 1-e^{-K^c \times T_s}$    ---> (is the minus sign right?)

I am looking for some conversion formula that can handle all cases (including single/multiple poles or zeros).

Note: I do not want to use any Matlab built-in functions like c2d(), zp2tf(), or tf().

there are at least 3, likely more, methods.

$z = e^{sT}$
impulse invariant
bilinear transform


Designing a sound that localises well I am trying to create or find a sound that localises well.

I read that emergency vehicles are changing their sirens; people found it difficult to localise the sound of the old siren (a sinusoidal sweep), and apparently this led to a lot of accidents. Apparently a single frequency is difficult to localise.

I've experimented using harmonic tones — I found that my brain can resolve higher frequencies, say > 600Hz, much better than lower frequencies,  say < 150Hz.

Unfortunately this is no good to me — I'm building a musical trainer game, and it's important that I communicate a direction without suggesting a pitch, so I need something else, maybe something like the sound of glass breaking or a clap -- anything that doesn't suggest a pitch.

I've experimented with static bursts, and found they don't localise nearly as well as the tones.

Could anyone suggest a good sample sound source?

I wonder whether something like a twig snapping would trigger some sort of primal reflex :)

You want something that has lots of features in the spectral domain and in the time domain. However you have to be careful what other content or meaning you would like to convey or not. Spectra conveys pitch and/or coloration and timbre. Time transients convey rhythms and/or patterns.

For example you could do band limited noise burst with different center frequencies and patterns. That would allow yo to vary localization, spectrum and rhythm without directly conveying pitch. You could also use different percussion samples with different rhythms.


fundamental question regarding bandwidth What I've learnt so far is that, if we change a time-domain signal into frequency domain, its bandwidth shows the range of frequency of radio waves that required to generate such signal.
However, for example in fiber optic,a 1550nm light is used to be the source and BPSK is used to modulate it, but if there is only a single frequency of light (1550nm) forming a signal, then what is the point of finding the bandwidth?

A single frequency is a pure sinusoid, and thus unmodulated (no BPSK).  Once you add any modulation, you no longer have a pure sinusoid.  The FT of this distorted sinusoid will have a non-zero bandwidth from the sidebands required to change the sinusoid from just pure to somewhat distorted (same as modulated).


Measuring length of DNA fibers from an image of single molecules I am a biologist with very little experience with image processing but have sufficient knowledge of MATLAB and have the image processing toolbox. Ideally I am looking for a MATLAB based solution, but an approach outlining how to go about it would also be helpful.

Update (28 Nov 2011)
It appears that there are certain problems (such as overlaps in signal and definition of color) when using composite images (which is what I presented in the initial question).  I am attaching separate images from the 2 channels: green  and red  (the turquoise regions in the composite image can be ignored), and the coposite image . The red channel is bad for 2 reasons: 1. It has poor contrast due to higher background, 2. Since Red seems to bleed into the green at the background level.

A feature is defined as a region on the composite image that has
Green-Red-turquoise-Red-Green or equivalently the 2 adjacent linear segments on the green and the red that are colinear and contagious.

I am hoping that looking at the images from the two separate channels makes identification of the features easier.

I have the following suggestions for the algorithm:

First identify co-linear green segments (and determine lengths of the green segments)
Determine if there are adjacent contagious and colinear segments facing toward each other (i.e. green->red->  <-red <-green) in the red channel. If yes define the red segment length from the point where green segments end  (because they will overlap with the green segments) till the point on the red segment that is closest to the other red segment of the feature. (i.e. one of the ends of the red segment is set to the end of the overlapping green segment).

Many thanks!

Background:

My question relates to extracting feature from an image:

The original image (tif) is located here:

  Image example 1 (dropbox)

This image is a composite of 3 channels (in tif format): red, green and turquoise. The turquoise colored fibers simply marks all the DNA we have on the coverslip. The feature of interest is the Green-Red--turquoise--Red-Green feature on the single DNA strand that is the middle of the image.

Red is generally the noisiest. This example is good because the contrast is good. However, sometimes the images are not so nice and there is hue throughout the image, so hard-coding a specific RGB value for the green and red color might not work for all images. Also, note that the fibers are not necessarily horizontal, they might be rotated (but never vertical).

Please see this image for an example:

The original image (tif) is located here:

  Image example 2 (dropbox)

Also, sometimes a single image has many such features and sometimes there are multiple features on the same DNA strand. Finally sometimes there might be only partial features (i.e. isolated green or isolated red or isolated green-red segments, but unpaired).

Question:

I would be grateful if someone can help me obtain the lengths of the individual segments of green and the red segments i.e. since the feature of interest is Green-Red--turquoise--Red-Green, each feature would have an array of 5 values (length of the first green segment, length of the first red segment, length of the turquoise segment, length of the second Red segment and the length of the second Green segment).

Example in Mathematica:

(* Get your image*)
img = Import["http://dl.dropbox.com/u/18072545/c_29.tif"];
(*Detect the extended minima and remove background*)
nB = ImageSubtract[img, ColorNegate@FillingTransform@ColorNegate[img]]

(*Separate RGB channels*)
cS = ImageAdjust /@ ColorSeparate[img]

(*Binarize*)
bcS = Binarize[#, .4] & /@ cS

(*Remove large elements*)
tH = TopHatTransform[#, DiskMatrix[2]] & /@ bcS

(*Detect lines using a Hough Transform*)
lines = ImageLines[#, .01, .8] & /@ tH
(*Plot them*)
Show[img, Graphics[{
   Thickness[.01], Red, Line /@ (lines[[1]]),
   Thickness[.006], Green, Line /@ (lines[[2]]),
   Thickness[.004], Blue, Line /@ (lines[[3]])}]]
(*Red and green are superimposed*)

Edit

Here you may see the Red and Green clusters apart. As you can imagine, you have to decide when a portion is red!


Exercise Related to Frequency Resolution and SNR I'm studying a book about DSP and trying to make exercises. Here is one I'm interested in:
A scientist acquires 65,536 samples from an experiment at a sampling rate
of 1 MHz.  He knows that the signal contains a sinusoid at 100 kHz.  He
needs to determine is if there is a second sinusoid at 103 kHz.  As a start,
he takes a 65,536 point DFT of the signal.  To his surprise, all he can see
in the spectrum is noise.  He estimates that the signal he is looking for is
4 times lower in amplitude than the random noise (i.e., SNR = 0.25). He also
estimates that the SNR will need to be at least 3.0 for the signal to be
detected, if present.  To improve the SNR, he decides to break the signal
into segments, and average their spectra.  Arrange your answers to the
following questions in a table.
a. If he uses a segment length of 16,384 samples, what will be the frequency
resolution (i.e., the spacing between data points) in the averaged spectrum?
Give your answer in hertz.  How many segments will be averaged?  What is the
SNR of the averaged spectrum?  Does this have the required frequency
resolution?  Does this have the required SNR?
b. Repeat for segment lengths of: 8192, 4096, 2048, 1024, 512, 256, and 128.
Here is my solutions:
1)
16384 samples
4 segments will be averaged
freq. resolution would be (25*10^4)/8193Hz ~ 30Hz
(it is told in the book that random noise reduces in proportion to the square-root of the number of segments)
snr = 0.25 *$\sqrt{4}$ =0.5 (not enough)
freq. resolution is enough to resolve between 100kHz and 103kHz sinusoid
2)
2048 samples
32 segments will be averaged
$10^6$/32 = 31250
freq. resolution would be 31250/1025Hz ~ 30Hz
snr = 0.25 *$\sqrt{32}$ = 1.4 (not enoght)
freq. resolution is ok
3)
128 samples
512 segments will be averaged
freq. resolution would be 1953/65Hz ~ 30Hz
snr = 0.25 *$\sqrt{512}$ = 5.6 (is ok)
freq. resolution is ok
Am I right or wrong? Where I can be wrong?
Thank you so much in advance!

Your calculation of the SNR is right and the DFT Bin Resolution is also OK.
One thing you're missing is the effective resolution due to "Windowing" effect.
The DFT of finite number of samples is interpolated by a Dirichlet Kernel (Like Sinc).  It means you resolution is also limited by the main lobe width of the Sinc which is proportional to the inverse of the window length.
Namely, using less samples means wider lobe which reduces the resolution and the ability to differ between two harmonic signals.


Signal filtering using Python I'm trying to create an application using python that is capable of recording an audio signal and detecting short glitches in the signal. For example, I could be recording a 1000Hz sine wave and looking for instances when the audio cuts out for a sample or two. I've been achieving this by removing all the 1000Hz frequencies using a notch filter and searching where the remaining peaks occur. My program seems to work fine when I begin recording a signal, finish recording and then filter the entire recorded signal, search for peaks etc but I would like to get it working in real time, where I can begin recording, analyze 5 seconds worth of audio then move on to the next 5 seconds of audio and so on until I finish recording. Unfortunately my notch filter that worked previously no longer seems to work correctly when I do this.
The code for my notch filter is shown below (I am using scipy for signal processing and pyaudio for recording audio))

f0 = 1000.0  # Frequency to be removed from signal (Hz)
Q = 30.0  # Quality factor
w0 = f0/(fs/2)  # Normalized Frequency
# Design notch filter
b, a = signal.iirnotch(w0, Q)
zi = signal.lfilter_zi(b, a)
z, _ = signal.lfilter(b, a, xn, zi=zi*xn[0])
z2, _ = signal.lfilter(b, a, z, zi=zi*z[0])
y = signal.filtfilt(b, a, xn)

The graph shown below is a filtered 5 seconds of audio. The large peak just after 1 seconds is a genuine glitch but the smaller peaks just after 2,3 and 4 seconds should not be there.

These smaller peaks also continue to occur when I analyze the next 5 seconds of audio, despite there being no glitches at all in this chunk of audio, shown below.

When I begin recording audio, finish recording and then check for glitches, these miniature peaks do not occur at all, they only occur when I try to divde the incoming audio into chunks and check for glitches as the audio is being recorded. Can anyone explain what they are (something to do with harmonics?) and whether or not I can get rid of them? Ideally I would also like to get rid of the peaks that occur at the beginning and the end of every 5 second chunk as well. I dont have much of a background in DSP and I'm still learning a lot about it so I would appreciate any help.

Hey I've been meaning to write this program, too.  :D

The way you use the initial conditions is to just pass them from one stage to the next.

Also lfilter returns (y, zf), not (zf, y).

So it should look more like:

zi = signal.lfilter_zi(b, a) * xn[0]
z, zi2 = signal.lfilter(b, a, xn, zi=zi)
z2, zi3 = signal.lfilter(b, a, z, zi=zi2)
z3, zi4 = signal.lfilter(b, a, z2, zi=zi3)

Though you don't need to rename each one, of course.  Initialize it with

zi = signal.lfilter_zi(b, a) * xn[0]

and then in each loop do something like

for x in chunks(signal):
    y, zi = signal.lfilter(b, a, x, zi=zi)
    look_for_glitches_in(y)


Low pass and High pass filter Coefficient Is there any way to tell whether a filter is high pass or low pass by observing only it's time domain samples or coefficients?

To elaborate a bit on Fat32's answer: the most straightforward thing to do is to compute (or estimate) the following two sums:

$$H(e^{j 0})=\sum_nh[n]\tag{1}$$

and

$$H(e^{j \pi})=\sum_n(-1)^nh[n]\tag{2}$$

where $(1)$ is the value of the frequency response at DC (i.e., $\omega=0$), and $(2)$ is the value of the frequency response at Nyquist (i.e., at $\omega=\pi$).

A low pass filter should have a relatively large value for $(1)$ and a very small value (ideally zero) for $(2)$. For a high pass filter the opposite is the case. If both values are small (and if $h[n]$ is not zero) then it's probably a band pass filter, and if both values are relatively large, it's probably a band stop filter. This of course only applies if you can assume that the filter approximates some standard frequency selective filter characteristic.


Why do we substract a background image and not divide it? Subtracting a background image for background correction is the standard way of doing background correction to obtain the foreground.

I am wondering if there are situations, where it would make sense to divide with the background image instead.

For example, it seems that having a uniform image intensity after division could be beneficial in interpreting an image. E.g in microscopy images.
It also seems physically more intuitive to have uniform baseline intensity after background correction/division than to have zero baseline intensity due to subtraction.

Division happens, indirectly, in cases where the image formation model is multiplicative. Concepts behind background varies on whether one only considers smooth fluctuations, or one includes noise. For the later, there is a lot of literature on multiplicative noise.

  important example(s]: speckle noise [...] proper shadows due to
  undulations on the surface of the imaged objects, shadows cast by
  complex objects like foliage and Venetian blinds, dark spots caused by
  dust in the lens or image sensor, and variations in the gain of
  individual elements of the image sensor array.

For a combined case, the principle of homomorphic filtering considers a nonlinear transformation on the image $f$ to cope with two main components (one slow, illumination $i$ and one fast,  reflection $r$). In the most basic case, $f=i\times r$ is passed through a logarithm before some linear filtering (and an inverse exponential on the results).This allows to reduce the noise without compressing the dynamic range too much.

It is suitable for instance to shading effects on surfaces, where one end is darker than the other, for instance, from Homomorphic filtering – part 2 (Mathworks):

In the mathematical morphology community, you can also find the LIP concept, for Logarithm Image Processing.

Multiplicative Noise Removal Using Variable Splitting and Constrained Optimization
Logarithmic Image Processing: Additive Contrast, Multiplicative Contrast, and Associated Metrics
Image restoration with multiplicative noise
Side work: [Multiplicative Calculus in Biomedical Image Analysis](https://doi.org/10.1007/s10851-011-0275-1


How mic setting will affect audio classification model? I have audio samples with three classes :

100 audio samples : class 'A'
100 audio samples : class 'B'
100 audio samples : class 'C'

Class 'A' and Class 'B' audio samples are recorded from one phone mic with same setting ( distance, volume etc ) but class 'C' samples are from different phone mic, ( each audio sample is from different phone mic )

I am working on a ML classifier to classify all three audio classes. My question is if I downsample all the audio signals to one frequency (i.e 16 kHz)

Will model still be biased because of different phone mic?
What affects model will face because of different phone mic and what is
the better solution for this problem other than downsampling to same frequency ?

We can't tell you what your classifier does, sorry.

But yes, you have a systematic bias in your data, and your classifier will cling to whatever is the strongest discriminator if it works as hoped. I will rename your classes to make this clearer:

100 audio samples : class 'Microphone 1 subclass 1'
100 audio samples : class 'Microphone 1 subclass 2'
100 audio samples : class 'Microphone 2'

If that strongest discriminator is the microphone used in the downsampled signals (and that is not an unreasonable assumption to make), then yes, your classifier performance will depend on the microphone used.

The "properest" way would be to go out and acquire some of class 'A' and 'B' with microphone 2; but assuming that is unreasonably much work: you could try playing back all sounds and record them through both microphones, and use random batches so that the you don't only see class C with the properties of microphone 2, but also classes A and B through 2, and C through 1. Maybe that is sufficient to reduce the role of the microphone characteristics for classification.

But really, you have a class that is really identical with being recorded via microphone 2, and your choice to call that class "C" instead of "recordings done with microphone 2" is your choice alone!


Phase difference calculation between two wave files I am doing comparison between two wave files, How to calculate phase difference between those files?
Second wave is generated by making some changes in the first one .I have tried calculating phase like this:-

y, sr = librosa.load('first.wav')
D = librosa.stft(y)
magnitude, phase1 = librosa.magphase(D)
y1, sr1 = librosa.load('second.wav')
D1 = librosa.stft(y1)
magnitude1, phase2 = librosa.magphase(D1)

than

calculated the difference(phase1-phase2).

Is this the correct way?

Is this the correct way?

This may be technically correct but won't give you any useful information. The reasons for this are actually quite simple. The phase between two signals can be computed
accurately only if your signals contain pure sine tones. You basically have only 2 scenarios where computing the phase difference between 2 signals makes sense.

a) the 2 signals both contain the same frequency component but have
a phase shift between them
b) the 2 signals contain each a
different frequency component

For case #a, the phase difference will remain the same throughout the entire length of the signals, whereas in #b, the phase difference will change as you traverse the length of the signals (it's assumed that the frequency in signal #2 is not a multiple of the frequency in signal #1 and the signals don't start in phase)

When you have a signal that contains multiple frequency components, the phase calculated by an FFT routine (including librosa.stft, which computes dft on multiple signal frames) is pretty
much random (you can think of it as being the average of all frequency components contributing to a particular frequency bin, including noise i.e. undesirable components).

This means that computing the phase difference for the  purpose of  finding the true phase difference between  2 random signals is not possible and is completely meaningless.
Computing the phase difference (within a single signal) can still be useful in other applications where it can be used to  more accurately estimate the "true" bin frequencies but here it's the relative phases that we're dealing with so it's different.

In conclusion, if your signals don't fit either case a# or case b#, you won't be able to do what you want and it doesn't even make sense to try to do that in the first place.


I have the following accelerometer signal, I want to segment or slice the signal into Red parts and Green parts I want to count steps which are basically represented by peaks in the signal. the red part in the signal doesnt correspond to any walking activity, I want to ignore them. Any help on segmenting or croping out the green parts of the signal.

Since I don't have access to your data, I found similar data online at https://archive.ics.uci.edu/ml/datasets/User+Identification+From+Walking+Activity#. It gives X-Y-Z accelerator readings. It is pretty clear that frequency content in the green and red sections are different. The green has significantly more high frequency components than the red. The frequency content changing over time had me think that short time Fourier analysis would be helpful.

So what I did was download the dataset and plot one of the signals to make sure we are looking at similar data. Here is a signal from the dataset I found online (pretty similar to yours!):

The next step would be to take a short time Fourier transform. MATLAB has already implemented this in the stft function. Using the defaults and plotting the result, I get this:

Now you can see, over time, how the frequency content of the signal changes. Just as we suspected, most of the signal has more energy in the higher frequency, and there are two spots where the high frequency content vanishes. Essentially, the problem becomes detecting the dark blue vertical lines in the image, or equivalently detecting the higher frequency components. Hope this helps!


How does one describe the “frequency response” of a mixer? Suppose one uses a vector network analyzer to measure the $S_{21}$ parameter of an RF mixer over a frequency range. Would it be appropriate to call this the “frequency response” of the mixer? I find this confusing because mixers are not LTI systems, and it is therefore meaningless to speak of the frequency response of a mixer.

But, from what I understand, experimentally measuring the $S_{21}$ parameter through a LTI system such as a lowpass filter, amplifier, or cable precisely yields the frequency response of the system. In that case, what does the $S_{21}$ parameter through a mixer (not LTI) correspond to?

A mixer is time-varying, but an ideal mixer is linear.  In radio usage, a mixer is periodically varying, which actually does give you an avenue to evaluate it in terms of a frequency response.
In the time domain, in a superheterodyne or direct conversion receiver, an ideal mixer takes a received signal $r_{rf}(t)$ and multiplies it by a sine wave to produce an intermediate-frequency signal $r_{if}(t) = A \cos \left (\omega_{lo} t \right)r_{rf}(t)$.  Note that the fact that this is a simple multiplication of the input by some independent quantity makes the system linear; the fact that the independent quantity isn't constant makes it time-varying.
You can take the Fourier transform of the above equation:$$R_{if}(\omega) = \mathcal{F}\left \lbrace r_{if}(t) \right \rbrace = \mathcal{F}\left \lbrace A \cos \left (\omega_{lo} t \right)r_{rf}(t) \right \rbrace$$
The Fourier transform of a multiplication is the convolution of the Fourier transforms of the multiplicands, so $$R_{if}(\omega) = \mathcal{F}\left \lbrace A \cos \left (\omega_{lo} t \right) \right \rbrace * R_{rf}(\omega)$$
The Fourier transform of a cosine wave at $\omega_{lo}$ is a pair of Dirac impulses at $-\omega_{lo}$ and $\omega_{lo}$, so the output of the mixer is, in turn, $$R_{if}(\omega) = \frac{A}{2} R_{rf}(\omega + \omega_{lo}) + \frac{A}{2} R_{rf}(\omega - \omega_{lo})$$
You can extend this for a mixer that multiplies by something other than a sine wave, such as a chirp or (commonly for actual RF circuits) a square wave.
So all of this was to get you ready for -- tada! -- the frequency response of a mixer.  The ideal mixer responds as above.  A real mixer has a frequency response that depends on what it does to the RF signal before it hits the mixing circuit, and what it does to the LO signal after it hits the mixing circuit.  So -- at this level of analysis -- you can expect the output of a more real mixer to actually look something like $$R_{if}(\omega) = H_h(\omega) R_{rf}(\omega + \omega_{lo}) + H_l(\omega) R_{rf}(\omega - \omega_{lo})$$
In the nomenclature I'm using within this answer, the $H_h(\omega)$ or  $H_l(\omega)$ would be your "frequency response of the mixer", depending on whether you were interested in the high-side or low-side mixing product.
The one thing that you can't do on a simple spectrum analyzer is just get the $S_{21}$ figure, at least assuming that the analyzer is looking at the inputs and outputs around the same frequency.  You need to input one frequency, and then look for the output at the correct shifted frequency.


Synthetic sounds to describe motion in 3d or 100d Watch a butterfly flitting about, or an optimizer chugging along in 3d or 100d --
a sequence of points $X_1 \ X_2 \ X_3\  \ldots$

How could one generate synthetic sounds that convey

moving fast or slowly
moving steadily in one direction, or oscillating ?

The question has two parts (which should perhaps be asked separately):

how to filter a high-dimensional path to a 1d signal that conveys the kind of motion
how to convert that to pleasant sounds.

The goal would be to listen to sounds made by a running optimizer
as it moves along in 3d or 100d,
and glean some info on how it's moving.

You say your points are strictly ordered:

  a sequence of points $X_1 \ X_2 \ X_3\  \ldots\ X_N$

Well, in that case, you can

either define "speed" as "euclidean distance (or any other metric) per in between two consecutive points":
$$s_n = ||X_{n+1} - X_n||,\quad n=\{1,\ldots,N-1\}$$
and interpolate a steady function $s(t)$ from that sequence, or
you can define speed as the first derivate of an interpolation function through all $X_i$. However, finding an interpolator for $N$ points in 100 dimensions isn't easy, and also, its utility would be pretty questionable, unless you have some goal in mind that clearly says why a specific interpolator is appropriate for your problem.

After you have a scalar value along a single axis, well, map to tone or loudness.

I really doubt the "illustrating" effect is overly great here: the choice of your interpolation will probably have the most important effect on what listeners will take as an impression from this, and that choice is not backed by the data itself – it's a model choice you make. Which would be ok, if you know something about the data (i.e. all the individual dimensions just contain harmonic oscillations would suggest classical low pass filtering would be an appropriate method to interpolate "speed" points), but if you don't and want to use this "audiolisation" to understand your data better, you'd end up fooling yourself. You will listen to your own assumptions, not the data.

However, you mentioned PCA, which indicates the points aren't actually to be understood sequentially! That of course makes things a lot harder. You'd essentially be looking for a $100-1$ hypersurface that describes how your points are mainly distributed in your $100$ dimensional space, and then look for the main gradient of that – and that leads to all the nice problems that you get with higher dimensional analysis (lack of "one" derivate of function, mapping values from something that isn't a vector space to a vector space, which introduces the need for maps on manifolds... meh).


Frequency and amplitude estimation with missing samples My input signal is given by $$ x[n] = A\sin(2\pi f n/F_s + \varphi)$$
The signal is uniformly sampled at a sampling rate of $F_s $[Hz]. However some of the samples are missing. Assume that there are not many missing samples and that their location is known (i.e. I know the timestamp for each sample),  I am interested in estimating $A$ and $f$.

My initial approach is to interpolate the vector and use standard methods (e.g. FFT). Another thought is to perform some kind non-linear regression (which I am unfamiliar with).

Any help would be appreciated.

I think best method for your is to use non-linear regression. These methods do not need to use equidistant data. I answered on similar questions (fft-bin-interpolation). Below is citation from this answer:

If your signal is quite similar to sinus with noise, then I think you can estimate frequency by using fitting or regression methods. Search Internet with keywords: "fit", "regression", "least square" and "sine". You can find a lot of program code, tools and so on. For example (I do not try it) - sine-function-fit See also good discussion in StackExchange "Given a data set, how do you do a sinusoidal regression on paper? What are the equations, algorithms?" There are a lot of good books and manuals for fitting/regression. I can recommend classical book  Numerical Recipes. See chapter "Modeling of Data". Especially  "Fitting Data to a Straight Line" for base understanding and  "Nonlinear Models" for understanding method of solution your problem. It is possible that there are specific methods for sine function.
Implementation of regression in your case can be effective. You must solve minimization problem for system by size of 3 (for amplitude, frequency and phase). Calculation of coefficients of this system is sums of some functions from all your points - O(N).


Does modified Morlet wavelet function satisfy the admissibility condition? Does modified Morlet wavelet function satisfy the admissibility condition ?
Do you have any reference for the answer please.
Thanks.

No, it doesn't (at least not in theory). You know that the admissibility condition is

$$\int_{-\infty}^{\infty}\frac{|\Psi(\omega)|^2}{|\omega|}d\omega<\infty$$

where $\Psi(\omega)$ is the Fourier transform of the mother wavelet $\psi(t)$. If $\Psi(\omega)$ decays sufficiently fast (which is the case for the Morlet wavelet) then the admissibility condition reduces to $\Psi(0)=0$, i.e. $\psi(t)$ must satisfy a zero-mean condition. However, this condition is not satisfied for the Morlet wavelet. So the Morlet wavelet is not admissible, but in practice this is no problem because $\Psi(0)$ is very small.

[1] M. Vetterli, J. Kovačević: Wavelets and Subband Coding


How can I track detected objects from frame to frame I want to track people's faces in a video stream (like a security camera), but I am running into a problem. I was using a Haar classifier to find the faces and sending the detected faces to a system to identify the person that was detected (so I can display their name under there face). The problem is that identifying a person takes more than a second, too slow and unnecessary to do every frame as the same people tend to stick around.

What is a good way to track people frame to frame?
I was thinking of something like storing every person, and in the next frame trying to find faces detected very close to where they were, but this does not work very well as faces can skip frames, move more than I expect, or "jump" faces swapping who is who.
Is there an accepted way to solve this problem?
What are some ways to overcome this?

Object detection is relatively a heavy task as you've notice. Detecting the object (in your case human face) in every and each frame would be cumbersome and computationally immense. Therefore, you need to employ an object tracking technique. There are various tracking algorithms, of which, KLT and mean-shift are the two popular ones.
KLT works based on object features (corner,..) but mean-shift (usually) is implemented based on intensity of image.
Personally, I recommend mean-shift for your application, because I think face color would be more unique in a image sequence than general features (corner,edge,..) of the face.
The only thing you might need to add, is verification of the tracker by your object detection algorithm, every like 50 frames to ensure sound functionality of the tracker.

For implementations see:

http://opencv-python-tutroals.readthedocs.io/en/latest/py_tutorials/py_video/py_meanshift/py_meanshift.html

http://www.pyimagesearch.com/2015/09/14/ball-tracking-with-opencv/


Ripple in time-domain using lanczos resampling I'm seeing ripple in the time domain when resampling a "unity signal", that is, all ones, with lanczos filters. The ripple/error exactly resembles a cosine with a period of the sampling rate - ie. the error is largest at a fractionate sample offset of 0.5. Using larger $a$ for the lanczos kernel reduces the error, but it never ceases to exist.

Is this supposed to happen, or am I making an implementation error somewhere? Doing the convolution by hand shows that the terms never sum up to 1 (unless at offset = 0). Truncated sinc shows the same problem.
Are there any windowed sinc filters that do not exhibit this kind of error? Otherwise, any other suggestions for resampling techniques? I'm generally looking for something that offers the performance of a lanczos kernel with $a$ > 10.

Edit: Thanks for the comments. I'll elaborate a little on what I want: So I'm interpolating a waveform display. What's important and interesting is, that the interpolation passes through the original sample points and resembles the 'true' shape of the signal as much as possible. This means this is good:

So under- and overshoot is not a problem in itself (technically, it's the whole point). This, however, is (note this is just a DC signal) - ie. nothing in the signal provokes the ringing, it is just an artifact depending on your fractionate interpolation offset:

Figure 1. Illustration of the problem with Lanczos2 interpolation: The sum of the four time-shifted Lanczos-window multiplied sinc functions (red) do not sum (black) to unity in the middle between horizontal coordinates 1 and 2.

A simple fix would be to distribute one fourth of the deviation from unity as function of the fractional position in the middle interval, with negative sign, to each quadrant of the Lanczos2 kernel. The equations are lengthy, but this is all very easily done numerically. The difference to the original Lanczos2 kernel is small:

Figure 2. Lanczos2 kernel (blue) and fixed Lanczos2 kernel (red).

Figure 3. Fixed Lanczos2 kernels (red) sum to unity (black).

The magnitude frequency response from the fixed Lanczos2 kernel looks just fine compared to the original, by some measures even better:

Figure 4. Magnitude frequency response of Lanczos2 interpolation (blue) and the same with the fixed kernel (red).

A similar fix could be tried on larger Lanczos kernels. I tested it on Lanczos3 and saw no problems.


Neural networks in system identification - What type of activation functions? I made a free software for all operative systems, even Android. It's called Deeplearning2C. It can train a neural network and generate C code and MATLAB-code. C-code for embedded systems and MATLAB-code for simulation.

I tried a example when I trained a neural network with linear data with constraints = nonlinear model in other words. It works very well. Here is the result when I use identity as activation functions.

I identification the model from the features (inputs)

$$-y(t-1), -y(t-2), ... , -y(t-n), u(t), u(t-1), ... , u(t-m)$$

And the label (output) $y(t)$. This is the regular transfer function or ordinary differential equation identification method.

But the issue I have is that I have a lack of experience of system identification when it comes to neural networks. I'm used to state space models (Subspace identification) and regular parameter estimation (Recursive Least Square). I wrote a library for that too.

My questions for you are that:

What type of activation functions should I use?
How many layers is "good" enough?
How many neurons should I have?
How do I know...how to know what to select? I'm seeking practical experience for identification with neural networks.

Here is my software Deeplearning2C. It's using Deeplearning4J as the core, but the user interface is made by me.

Also here is the MATLAB/GNU Octave library for linear system identification with recursive least square and subspace identification. I wrote a C-library too for recursive least square identification as well.

Looks like you've done a lot of work on your projects. As @MarcusMüller said, by far the majority of people start with ReLU and go from there. It doesn't have the "vanishing gradient" problem that tanh has for example. All your questions are open ended but common for designing neural networks. There are so many "nobs to turn" to try and make your network be the best it can be.

You could adjust some parameters, train, test, and compare. But then days or weeks have gone by and you'll still be tuning your parameters. One popular method to choosing "good" parameters is called hyper parameter optimization. This blog post does a good job explaining it: https://blog.floydhub.com/guide-to-hyperparameters-search-for-deep-learning-models/. The idea is that you randomly sample your parameters space (for example, number of neurons, a few different activation functions, different learning rates, number of epochs, etc.) and in search for the a good set of parameters. I hope this helps!


DFT of 2d real signal and Hermitian symmetry Knowing that DFT of n-values real signal in 1d consists of n/2+1 different values where the second half of the spectrum is complex conjugate of the first one (Hermitian symmetry).

However in the spectrum of 2d signal that I've made in Octave I see complex conjugates in both dimensions. However the FFTW manual states that we can take a use of the symmetry in only one dimension.

What is the truth? Is Hermitian symmetry only useful in one of multidimensional transform or could it be used in other dimensions?

[EDIT]

Ok, so since the Hermitian symmetry occurs in each of the directions then why the transformed vector in FFTW is diminished only by half (only one dimension values are stored on N/2+1 complex numbers)? Why isn't each direction diminished?

Yes, multidimensional DFT's have complex conjugate symmetry in each dimension

I think you might be misinterpreting the FFTW manual which is referring to internal storage formats.

Edit:

  why the transformed vector in FFTW is diminished only by half (only
  one dimension values are stored on N/2+1 complex numbers)? Why isn't
  each direction diminished?

From the FFTW manual:

  The multi-dimensional transforms of FFTW [...] compute simply the separable product of the given 1d transform along each dimension of the array

So for example if you provide a 256x256 array of real numbers the result will be 256 1 dimensional DFT's. As they mention at the bottom of that documentation this is not a multi-dimensional DFT.

see http://www.fftw.org/fftw3_doc/Complex-Multi_002dDimensional-DFTs.html


Probability of False Alarm versus SNR I'm currently testing the performance of an energy detector using two tests. One test involves finding the optimal threshold to achieve a probability of detection(Pd) of 0.9 and a probability of false Alarm (Pfa) of 0.1.

Having completed the first test, the next test I want to do is to see how varying SNR affects both Pd and Pfa, given a constant threshold. If Pd is the probability of the detector detecting a signal given the signal of interest is being transmitted and pfa is the probability of the detector detecting a signal given the signal of interest is not being transmitted how exactly can I vary SNR to get Pfa. I'm sure its something that I'm not quite understanding but if the signal of interest is not being transmitted and all that the receiver is picking up is noise, how is there a SNR in the case of Pfa?

Hi: In statistics we call that the probability of type I error ( rejecting when true ) and type II error ( accepting when false ). The way it's done there is that, once you make an assumption about the distribution of the data ( i.e: normal, t, whatever ), you decide on the null and the alternative, along with what you want the P(type I error ) to be ( say 0.05 ) and then, given the assumption about the distribution, your type II error probability is then known or can be calculated ( sometimes with difficulty ). Your case may be slightly different but I would think that a statistical signal processing book would have a "hypothesis testing" section which would explain what I am describing in more detail.  I hope this helps a little. I put it in the answer section even though it's not really an answer because it was too long for a comment.


Determine proper coefficients for the transfer function The problem says that we have the information (actually a .mat file that contains the data) of 1 pixel from an IR camera from an object with 26° of temperature. Because of the noise, the output is given by

$$
y[n]=26+N[n]
$$
where $N[n]$ represents the noise. For that, I'm supposed to implement the following filter:
$$
H[z]=\frac{b}{1-az^{-1}}
$$
where $ a\in (0,1)$ and $b\in \mathbb{R}$. The directions say "fix a value for $a$ then find a proper value for $b$ such that you get the known value of temperature". I'm working on Matlab and it seems the best values are $a=b=0.5$ (just by inspection by looking the graphs after the filter). How can I exactly determine a good value for $b$?

Your are basically designing a lowpass filter in order to suppress the noise to some extend. For now I will assume that you are dealing with white noise, so a flat frequency spectrum with a time average of zero. But the value you are after is 26, which in the frequency spectrum is just a spike at zero Hertz. So in order for your filter to recover this value a good as possible it should be stable and have unity gain at zero Hertz.

Fixing $a\in(0,1)$ should ensure that the transfer function has a pole inside the unit disk and therefore be stable. In order to ensure the unity gain at zero Hertz you need to evaluate $H(z)$ at zero Hertz. In order to evaluate a discrete time transfer function at a frequency you can use $z=e^{j\,\omega\,T}$, with $\omega$ the frequency is radians per second and $T$ the sample time. So for zero Hertz/radians per second you get $z=1$. So evaluating $H(z)$ at this $z$ value gives

$$
H(1) = \frac{b}{1 - a}.
$$

This is only equal to one (unity) if $b=1-a$. So if you used $a=0.5$ then indeed it should be that $b=0.5$.


How to detect the edges of periodic signal? I have a periodic signal generated by accelerometer sensor and I need to count the number of peaks. However, I can't predict when the periodic sequence starts and ends. What would be the best way to do this if expected frequency of the useful signal is between 0.2Hz - 5Hz?

Here is the example of the signal - five squats.

The autocorrelation function of a sliding window of data may change by a significant amount when the window is slid past the first or last cycle of your burst of periodic events.  You may have to play with the width of the autocorrelation window depending on the estimated frequency of your cyclic/periodic events.


Machine Learning in FAST corner detector This paper by Edward Rosten and Tom Drummond 2008 describes choosing interest point via machine learning approach using ID3 algorithm. I understand the concept of the ID3 they are using. But my question is how are they fitting this whole training and testing idea of machine learning, while are analysing real time video frame.

But my question is how are they fitting this whole training and testing idea of machine learning, while are analysing real time video frame.

Your question is not clear, but I'll try to answer. They don`t analyze "real time video frame" with machine learning. What do they do:

They have slow segmented test == test all points around whole circle, if it looks like corner. No machine learning.
They want to do it fast, so to test not all points around circle, but only some of them. And they want to know which points test first to have answer.
That`s why they take some training data and use "slow" test in it - to have ground truth. Then they train decision tree to approximate this "slow" test.
Basically - you have (for simplicity) 8 points for test around circle. Each point can be darker, brighter or the same as central. Result of testing can be "corner" or "not a corner".
Example. You test north point, it is "darker". You look for ground truth and see, that there are both "corner" and "not a corner" points with "darker" north point. Cannot decide. Then you test south point - it is also darker, then east and west - all are darker.You look for the ground truth and see that all points with this test results are "not a corner".
Now you can answer "not a corner" with only 4 tests instead of all 8. So, you can save it as one branch of decision tree. The same for other cases. You find optimal way by some criterion and form full decision tree.
You or me use FAST detector - we don`t retrain trees, we just use ready ones to have it, so no machine learning in real time.


FSK carrier frequency estimation In the MPSK Modulations,  Cyclic Spectrum Density method can give a peaks, for example for the QPSK:

But in the FSK signal, this peaks do'nt appear. How to solve it?

To perform carrier recovery for FSK, if the data is random and can assumed to be equiprobable, then you can take the mean of the derivative of the phase versus time to determine the carrier offset.

Optionally the phase change from one sample to the next can be estimated using using a complex conjugate product of successive samples (for small angles, the imaginary result of the complex conjugate product of two normalized vectors is approximately the phase). Since this is the change in phase over a time interval, this result will be directly proportional to frequency:

$$f \propto Im(s(n)s(n-1)^*)$$

Where s(n) is the normalized complex FSK time domain signal.

Note in terms of I and Q samples this is simply:

$$I_{n-1}Q_n-I_nQ_{n-1}$$

The average of this result is proportional to the frequency of
the carrier and the deviation from average is the demodulated FSK signal. (The carrier component assuming a fixed offset can therefore be removed by using a simple DC Nulling Filter such as the one shown in this post: What does correcting IQ do?)

Since a small angle approximation was used, this provides a linear demodulation when the sampling rate is sufficiently high such that the phase change between two samples is a small angle such that $sin(\theta) \approx \theta$. For FSK that linearity is generally not as critical so would also work with much larger phase angles.


Downsample: resample vs antialias fitlering + decimation I have a discrete signal sampled @Fs. I need to downsample it to Fs/k. Main reason is to reduce signal bandwidth and speed up computation : I'm only interested in a reduced bandwidth < Fs/2k. Some high frequency noise >> Fs/2k can exist : antialiasing is needed.
From what I understand, I have several options :

brute force : decimation (without taking care of aliasing problem!)
antialiasing + decimation : I can use for example numpy
decimate function
frequency approach : FFT then truncation then IFFT, can be done
    using resample

Only solution 2- 3- help solving antialiasing issue.
3- filters in frequency domain, supposing signal is periodic. 2- filters with IIR/FIR filter.

I can always use decimate without questioning, but I would like to precisely understand limitations behind. What are drawbacks / advantages of using solution 2- versus 3-?
Do you have some practical guidelines about downsampling? Some recommendations?

I will explain why method 2 is often a better choice over method 3.
The frequency domain approach is equivalent to the "Windowing" method of filter design- in that to do that approach correctly you should window your data before taking the FFT. For an anti-alias filter design in the time domain approach, the least squares filter design algorithm outperforms window design approaches. (See this post for a detailed discussion on that: FIR Filter Design: Window vs Parks McClellan and Least Squares). For time-domain filters for decimation and interpolation applications, the least-squares filter design is a better choice over equiripple because of the stop-band roll-off: for equiripple the stop-band is at the same level in each aliased frequency band resulting in more overall noise folding in than you would get with least-squares.
Side Note: If you observe the coefficients for an equiripple design, you will often observe (if the filter isn’t too long) two slightly larger "impulses" toward the beginning and end of the filter impulse response (the coefficients of the filter is the impulse response). Remove those larger coefficients near the tails of the response and the equiripple design will also have the desired feature of stop-band roll-off! For further details on that please see: Convert a Park McClellan FIR Solution to Achieve Stop-band Roll-off
Further the least-squares (and equiripple) design tools in Matlab/Octave/Python feature multi-band filter design which is ideal for decimation (and interpolation) applications since the images are limited to distinct bands. Thus you can optimize the filter rejection to just the frequency locations that would fold in, further optimizing the solution given the same number of taps. Below is an example spectrum I have recently shown for interpolation and the resulting multi-band filter designs for both least squares and equiripple appropriate for eliminating the images (this is the interpolation filter to grow the zeros that are inserted to their interpolated value by eliminating the images, the same would apply to the decimation filter where we want to reject these same image locations prior to throwing samples away). This also converts readily to an efficient polyphase filter structure by mapping the same coefficients row to column in the polyphase filter, as detailed in other posts here.
In this plot the blue is the desired spectrum along with its images, and the red and the black show the multi-band filter response for the two different filter design choices (red is least squares and black is the Parks-McCelllan or equiripple design). This would be equivalent for a decimator except the images would be noise or other signals that could fold in during the process of throwing away samples for decimation. Given the same number of taps observe how the total noise that would fold into band is significantly less with the least squares filter design.

Note: If you have enough samples so as to not truncate the desired response, you could certainly still do the least squares filter design approach in the frequency domain ---- the filtering (convolution) described above that is done in the time domain is equivalent to multiplying in the frequency domain- but to do this properly would necessitate a lot more samples to ensure sufficient tails of the kernel (the frequency transform of the filter's impulse response) are included.


FIR filter phase response With reference to the diagram below, is it correct to say that at the pass band the phase response is linear. It also remains linear in the stop band but it get chunks of π. These jumps happen when the magnitude response goes down to minus infinity dB (theoretically) ?

You are right. Even though the phase jumps at the zeros of the frequency response, such a phase response is usually still called "linear". For a frequency selective filter with frequency response zeros in the stopband, the phase always has discontinuities at the zeros. A purely linear phase response (without jumps) is only possible for filters with no zeros in their frequency response.

The frequency response of your filter can be written as

$$H(e^{j\omega})=|H(e^{j\omega})|e^{j\phi(\omega)}\tag{1}$$

where $\phi(\omega)$ is the phase response with jumps at the zeros, as shown in your figure. For such a linear phase system, the frequency response (1) can equivalently be written as

$$H(e^{j\omega})=A(e^{j\omega})e^{j\hat{\phi}(\omega)}\tag{2}$$

where $A(e^{j\omega})$ is a real-valued but bipolar (i.e. positive and possibly negative) function satisfying $|A(e^{j\omega})|=|H(e^{j\omega})|$. The phase $\hat{\phi}(\omega)$ in (2) is now a purely linear function without any jumps.


Why does each node in the wavelet scattering transform split into multiple paths? Why does each node in the wavelet scattering transform split into multiple paths as in this figure from https://arxiv.org/abs/1304.6763?

I understand roughly what’s happening along a single path, but I don’t understand why each node splits into three.
Thanks!

We take a wavelet transform over other wavelet transforms, which involves multiple wavelets. E.g.:

First-order uses 100 wavelets to tile the frequency domain and produce 100 coefficients
Second-order uses 16 wavelets to tile the frequency domain and produces 8 coefficients per first-order coefficient, yielding a total of 800
Third-order uses 9 wavelets ... etc

(Note second-order wouldn't actually produce 800 in implementation since xi2 > xi1 wavelets are discarded, where xi == center frequency).
The idea is, first-order coefficients capture temporal modulations (AM-FM). Second-order capture the modulations of those modulations, extracting higher order structure, and the process repeats. Progressively we require less wavelets as there's less fine-grained structure to capture (also energy and non-redundancy reasons), but it's still more than one wavelet (i.e. path), as just one wavelet cannot cover entire frequency domain without losing information (unless it's very wide, but that's a bad wavelet).
See helpful lecture, and Python implementation. I also might write a detailed post on this topic in coming weeks.


Differences in moving window rms due to sampling rate I have a 10 second recorded signal band pass filtered from 100-500 Hz. The original sampling rate is 10 kS. Now I run a moving window RMS with a window length: 10 ms on the signal (10/1000 * sampling rate). I then find the samples of the RMS signal that are greater than 2 SD from mean. I get a total number of, say, 4894 samples that satisfy this condition.

Now I down sample the signal to 2 kS, and do the same process with the recalculated window length also of 10 ms (but different sample length due to the sampling rate).
When I apply the same >2 SD criteria I get 974 samples.

I then use the sample indices from the two different cases to see if they are indeed the same samples by normalizing the indices by sampling rate times 1000 to get ms precision.

Then I compute the intersection of these two sets of times when s.d > 2 for both the sampling rates.

I get a value of around 700. I am trying to understand why. Shouldn't all the samples times in the set corresponding to the 2 kS sampling rate actually be a subset of the sampling times seen in the 5 kS sampling rate for the > 2 SD criteria?

Shouldn't all the samples times in the set corresponding to the 2 kS
  sampling rate actually be a subset of the sampling times seen in the 5
  kS ?

Not necessarily. It depends on how you do the down sampling. If you use an anti-aliasing lowpass filter, it typically induces time delay and/of phase distortion. If you don't, you probably end up with some residual aliasing, since there is no such thing as a true "band limited" signal.

Could also be rounding errors when you determine the sample locations.


How do I manually plot the frequency response of a bandpass Butterworth filter in MATLAB without freqz function? I have code like below that applies a bandpass filter onto a signal. I am quite a noob at DSP and I want to understand what is going on behind the scenes before I proceed.

To do this, I want to know how to plot the frequency response of the filter without using freqz.

[b, a] = butter(order, [flo fhi]);
filtered_signal = filter(b, a, unfiltered_signal)

Given the outputs [b, a] how would I do this? This seems like it would be a simple task, but I'm having a hard time finding what I need in the documentation or online.

I'd also like to understand how to do this as quickly as possible, e.g. using an fft or other fast algorithm.

We know that in general transfer function of a filter is given by:

$$H(z)=\dfrac{\sum_{k=0}^{M}b_kz^{-k}}{\sum_{k=0}^{N}a_kz^{-k}} $$

Now substitute $z=e^{j\omega}$ to evaluate the transfer function on the unit circle:

$$H(e^{j\omega})=\dfrac{\sum_{k=0}^{M}b_ke^{-j\omega k}}{\sum_{k=0}^{N}a_ke^{-j\omega k}} $$

Thus this becomes only a problem of polynomial evaluation at a given $\omega$. Here are the steps:

Create a vector of angular frequencies $\omega = [0, \ldots,\pi]$ for the first half of spectrum (no need to go up to $2\pi$) and save it in w.
Pre-compute exponent $e^{-j\omega}$ at all of them and store it in variable ze.
Use the polyval function to calculate the values of numerator and denominator by calling: polyval(b, ze), divide them and store in H. Because we are interested in amplitude, then take the absolute value of the result.
Convert to dB scale by using: $H_{dB}=20\log_{10} H $ - in this case $1$ is the reference value.

Putting all of that in code:

%% Filter definition
a = [1 -0.5 -0.25]; % Some filter with lot's of static gain
b = [1 3 2];

%% My freqz calculation
N = 1024; % Number of points to evaluate at
upp = pi; % Evaluate only up to fs/2
% Create the vector of angular frequencies at one more point.
% After that remove the last element (Nyquist frequency)
w = linspace(0, pi, N+1);
w(end) = [];
ze = exp(-1j*w); % Pre-compute exponent
H = polyval(b, ze)./polyval(a, ze); % Evaluate transfer function and take the amplitude
Ha = abs(H);
Hdb  = 20*log10(Ha); % Convert to dB scale
wn   = w/pi;
% Plot and set axis limits
xlim = ([0 1]);
plot(wn, Hdb)
grid on

%% MATLAB freqz
figure
freqz(b,a)

Original output of freqz:

And the output of my script:

And quick comparison in linear scale - looks great!

[h_f, w_f] = freqz(b,a);
figure
xlim = ([0 1]);
plot(w, Ha) % mine
grid on
hold on
plot(w_f, abs(h_f), '--r') % MATLAB
legend({'my freqz','MATLAB freqz'})

Now you can rewrite it into some function and add few conditions to make it more useful.

Another way (previously proposed is more reliable) would be to use the fundamental property, that frequency response of a filter is a Fourier Transform of its impulse response:

$$H(\omega)=\mathcal{F}\{h(t)\} $$

Therefore you must feed into your system $\delta(t)$ signal, calculate the response of your filter and take the FFT of it:

d = [zeros(1,length(w_f)) 1 zeros(1,length(w_f)-1)];
h = filter(b, a, d);
HH = abs(fft(h));
HH = HH(1:length(w_f));

In comparison this will produce the following:


QPSK and OQPSK Modulation In QPSK, after splitting the odd and even bits into the I and Q channels, and then multiplying them with the sine and cosine function, what next am I suppose to do?

Do I rearrange them back into their odd and even positions?

Also, in OQPSK do I make the first bit of my Q array a dummy bit zero to represent the half symbol delay in the Q channel ?

This question seems to be based on several misconceptions.

In a QPSK modulator operating at a rate of $N$ baud, two bits enter the
modulator during each $T = N^{-1}$ second interval.  If the bits are
entering on a single wire (at a rate of $2N$ bits per second, the
QPSK modulator first converts this serial input
bit stream at $2N$ bits per second into two parallel bit streams at
$N$ bits per second. For example, alternate bits are steered to
alternate streams so that, say, all the odd-numbered bits
are in one stream and all the even-numbered bits are in the
other stream.

In some cases, the bits entering the QPSK modulator
are already formatted into two parallel bit streams at $N$ bits per
second. One important and commonly observed arrangement of this type
is when a rate-$\frac12$ convolutional encoder precedes the QPSK modulator,
so that information bits enter the encoder serially
at rate $N$ bits per second, and the convolutional encoder spits out two
parallel bit streams at $N$ bits per second.  However, regardless of
encoding or not, the bits that the QPSK
modulator uses are in two parallel streams that are synchronous
(that is, the bit transitions occur at the same times (once every
$T$ seconds) and in each
bit stream, the bit occupies the entire $T$ second interval.

As the OP's question says, the bits are modulated onto cosine and
sine carriers.  As discussed in more detail here, if the bits in the two streams are denoted by $b_I$ and
$b_Q$, then the corresponding signals, say, over the interval
$[0,T)$ can be taken to be
$$(-1)^{b_I}\cos(2\pi f_c t) ~~ \text{and}~~ -(-1)^{b_Q}\sin(2\pi f_c t),
0 \leq t < T. \tag{1}$$
(Ignore the negative sign on the $\sin$, it is there to simplify
notation in more detailed analyses).  The QPSK signal is the sum
$$(-1)^{b_I}\cos(2\pi f_c t) -(-1)^{b_Q}\sin(2\pi f_c t),
0 \leq t \leq T. \tag{2}$$
It is not the case, as the OP seems to think that the two
signals are separated into even and odd bits, so that
$(-1)^{b_I}\cos(2\pi f_c t)$ lasts from $t=0$ to $t=T/2$ and
$-(-1)^{b_Q}\sin(2\pi f_c t)$ from $t=T/2$ to $t=T$: the
cosine and sine signals last for the full $T$-second duration
and the bits $b_I$ and $b_Q$ both affect the signal phase for all
$T$ seconds.

Offset QPSK also converts the serial bit stream into
two parallel bit streams in which the bit intervals are
of duration $T$ on each wire, but the bit transitions
are not synchronous but occur with an offset of
$T/2$ seconds.  Thus, $b_I$ modulates the cosine carrier
to produce
$$(-1)^{b_I}\cos(2\pi f_c t), 0 \leq t < T$$
but $b_Q$ modulates the sine carrier to produce
$$-(-1)^{b_Q}\sin(2\pi f_c t), \frac{T}{2} \leq t < \frac{3T}{2}.$$
The modulator output is still the sum of the two sinusoids,
but the sinusoids change phase at different times. In fact,
the signal can be deemed to be a $2N$ baud signal since the
phase can change once every $T/2 = (2N)^{-1}$ seconds.
Each bit still affects the phase for $T$ seconds, but the
bit transitions are staggered.  In the receiver, integrations
(or matched filterings)
are carried out over (staggered) $T$-second intervals in the
two branches of the receiver, and the decision devices
spit out one bit every $T$ seconds but at staggered times.
It is easy to multiplex these decisions onto a single wire
that thus carries a bit stream at $2N$ bits per second,
just as the input to the QPSK modulator (if anyone still
remembers that far back in this answer) is a bit stream
on a single wire at $2N$ bits per second.


What is the "bilinear interpolation kernel" in personlab paper? please excuse my ignorance in computer graphics, but what is this bilinear interpolation kernel in the personlab paper page 6 equation (1)?

Here it is:
$$
h_k(x) = \frac{1}{\pi R^2}\sum_{i=1:N}p_k(x_i)B(x_i+S_k(x_i)-x)
$$
where:

$x$, $x_i$ are image coordinates (2-D vectors).
$h_k \colon \mathbb{R}^2 \mapsto \mathbb{R}$ localized heatmaps obtained by a sort of Hough voting.
$S_k \colon \mathbb{R}^2 \mapsto \mathbb{R}$ offest vectors
$p_k \colon \mathbb{R}^2 \mapsto \mathbb{R}^2$ heatmaps
$B(.)$ denotes the bilinear intepolation kernel.

  This is a
  form of Hough voting: each point i in the image crop grid
  casts a vote with its estimate for the position of every keypoint, with the vote being weighted by the probability that
  it is in the disk of influence of the corresponding keypoint

I guess bilinear interpolation kernel is a class of function and not an actual one (I may be wrong).

But if it is bilinear why does it only take one argument in the equation?

Maybe this is bilinear 1D interpolation, intepolating between the coordinates of 2-D vector?

Otherwise this formula doesn't make sense...

When I read «bilinear interpolation kernel» in image processing, I assume that they mean bilinear image imterpolation:
y(n+t) = (t-1)x(n) + tx(n+1)

If you decide upon a fixed uniform upsampling factor and standard dsp «zero-fill then lowpass filter» upsampling, the convolution kernel turns out to be samples of a triangular waveform.

-k


Signals and systems book with linear algebra approach I am currently taking a course on discrete-time signals and systems. I am using B. P. Lathi's Signal Processing and Linear Systems, which I don't like at all, as it doesn't draw any paralells to linear algebra. The concept of eigenvalues and eigenvectors in linear transformation is never discussed, for example.

Can you please recommend a book on signals and systems with a linear algebra approach?

I'm not particularly current, my junior level signals and systems text was

  Lathi, Bhagwal. Signal, Systems, and Controls. Intext, 1973

which had a some linear algebra, primarily because it had a control engineering bent.  The last chapter on mixed continuous-discrete systems was a lot to digest. I have come to appreciate it more over the years.

I used an older version of Chen in graduate school and it is very much from a linear algebra perspective, but is not aimed at the same audience Lathi is.

  Chen, Chi-Tsong. Linear system theory and design. Oxford University
  Press, Inc., 1998.

Kailath is a bit harder than Lathi,

  Kailath, Thomas. Linear systems. Vol. 156. Englewood Cliffs, NJ:
  Prentice-Hall, 1980.

and Luenberberger is very readable.

  Luenberger, David G. Introduction to dynamic systems: theory, models,
  and applications. Vol. 1. New York: Wiley, 1979.

In signal processing,

  Oppenheim, Alan V., and Ronald W. Schafer. Discrete-time signal processing: Pearson new International Edition. Pearson Higher Ed,
  2013.

Is sort of the "standard" text book on signal processing but light on linear algebra.

Strang is good but assumes you have some linear systems background.

  Strang, Gilbert, and Truong Nguyen. Wavelets and filter banks. SIAM,
  1996.

and his book on Linear Algebra is a very good introduction.

Its hard to recommend a book without knowing your back ground, but If I were to ask the same question you did,  I would go with Chen.


FFT processing of IQ data My time-domain sampling rate of IQ data is Fs. Matlab FFT assumes spectral bandwidth is Fs/2, per the Nyquist criterion. I understand that IQ data has bandwidth of Fs. How do I do FFT of such IQ data and get bandwidth of Fs?

For strictly real data (imaginary components all zero), if you only look in the first half of an FFT result (N/2 result bins for an FFT of length N) you only see Fs/2 bandwidth.  The second half of that FFT result is just a redundant conjugate mirror of the first half given strictly real input.

But for complex input or IQ sampled input, the 2nd half of the FFT result is not just a conjugate mirror of the first half, but represents an independant negative spectrum, of bandwidth Fs/2.  Total bandwidth of both halves of the complex input FFT being Fs/2 + Fs/2 = Fs (minus the roll-off or transition region of any anti-aliasing filter before IQ sampling).

For direct conversion from f0 to baseband IQ (say using a Tayloe mixer or other quadrature heterodyne) sampled at Fs (then doing an FFT on the IQ result), the spectrum from f0-Fs/2 to f0 is in FFT result bins N/2 to N-1 plus bin 0, and the spectrum f0 to f0+fs/2 is in FFT result bins 0 to N/2.

The Nyquist criterion is sort of met because an IQ sample vector really contains twice as many samples (real components plus quadrature components), thus providing twice the bandwidth. But a complex vector can use the same length of a complex FFT as a strictly real vector of the same time duration.  So conflating the FFT length with the number of real samples or complex samples or components of complex samples can be confusing.


Finding Distance between Cameras in Stereo Vision We know how to create disparity map from stereo vision with OpenCV. My question is How can we find/estimate distance between cameras from two image taken from these cameras in OpenCV?

Calibrate your cameras with the same world coordinate system. Use the same checkerboard and take images at the same time with both cameras.

After you have your extrinsic matrices for each camera, you can use this answer here where you need to solve the equation for the position of each camera center in world coordinate system.

Then you can calculate the distance between two coordinates.