Channel estimation and interpolation by vector in OFDM I have V-OFDM system (Vector OFDM) where every symbol is organized as in below figure:

The black cirles are pilots and white circles are data. First, iFFT operation it taken per column, and then collecting data is taken row-wise resulting collected data as shown in the figure right. It means first eight samples are pilot and following 24 sample are data, then 8 samples are pilot and so on. Then adding CP --> channel --> noise -->received data
In the receiver side, After usual processing of removing CP -->reshaping --> FFT column wise. We get the received data equivalent into transmitted data in above figure left. Then extract the pilots which are in our case of size $(3,8)$.
Then, LS estimation is performed, by dividing the received pilot by original pilot to get the frequency-domain estimated channel.
I ask here, how can we interpolate the channel in this case ? if I get a frequency-domain estimated channel of size $(3,8)$, what's the process to interpolate it into size of my symbol ? I think it should be vector interpolation? how can I do it in that case?
Thank you

You can use various methods to interpolate the channel - Linear, Polynomial, Sinc Interpolation etc. But what you need to keep in mind is synchronization. You have to make sure that frequency and timing offsets are eliminated or accounted for. Otherwise you will see an error floor in the Channel Estimation Error. Means your Mean Square Error (MSE) for Channel estimate will not go down even if SNR increases. It will go down for lower SNR but will flatten out for higher SNR. The residual errors due to timing/frequency offsets will spill over to channel estimation. See this paper (https://ir.nctu.edu.tw/bitstream/11536/32821/1/000072736100027.pdf). A simple search in google scholar with 'Channel Interpolation in OFDM' will reveal lot of good references.


Produce PSD from a FFT output I am using an FFT function in PHP language (https://gist.github.com/mbijon/1332348) which works very well, but cant produce PSD from its output:
$fft = Fourier($data,1);
// $fft[0] is real portion of 1
// $fft[1] is imaginary portion of 1
// $fft[2] is real portion of 2
// $fft[3] is imaginary portion of 2
// ...
// $fft[n-1] is real portion of n
// $fft[n] is imaginary portion of n

PSD formula In Matlab:
PSD = Y.*conj(Y)/N

I cant understand this formula in math notation form...
How can I adapt this formula to work with that output?

Indeed, PHP is not a weapon of choice to do signal processing. But if you absolutely need to "adapt" the matlab formula, here is an one-liner to "conjugate" the data, based on commented lines from the code of your question (// $fft[0] is real portion of 1, etc.):
// PHP's $conjdata[·] stands for matlab's conj[Y]
for ($i = 0, $cnt = count($data), $ri = 1; $i < $cnt; $i++, $ri = -$ri) $conjdata[$i] = $data[$i]*$ri;

Also, the PSD formula is somewhat trickier than just Y*conj(Y)/N. See, for example, https://stackoverflow.com/questions/3949324/calculate-autocorrelation-using-fft-in-matlab, the second answer.


Applying a Power on a Sum of Gaussians to Enhance Resolution There is so-called Power Transform technique in signal processing where the unresolved signal is sharpened by raising each data to a constant positive power. The example shown below is a sum of two Gaussians, $G_1$ and $G_2$ plotted against time. In order to resolve the peaks one can apply a power of 5 to each data point. The result is the resolved pair as shown in the bottom figure.

The above operating is equivalent to raising the sum of Gaussians $G_1$ and $G_2$ to the power 5, i.e., $(G_1 + G_2)^5$ then according to the binomial theorem, we should get 6 peaks, because not only we have the first and the second peak raised to power 5, we will have a "mixture" of the two peaks as well- the other terms of the binomial theorem. No matter how well we zoom, we cannot see the missing $G_1^4$$G_2$, $G_1^3$$G_2^2$ peaks? If power transform of 5 is applied to two peaks, we are left with two peaks not six.

Where do the peaks corresponding to the other binomial terms go?

Thanks.

The binomial theorem is not necessarily involved: the top waveform is simply pointwise raised to the 5-th power, as others have noted, and there are no missing peaks that should be seen. To illustrate, consider the following crude hand-drawn figure:

I drew this in my simulation software, so it was automatically "digitized" as I drew it. Obviously, I am poor at freehand drawing with a computer mouse, but I was only trying to make two peaks, with no underlying peak shape in mind.

Then the second figure shows what happens if the digitized values from the first figure are raised to the 5-th power and plotted:

Still ugly, but better resolved.

Now taking into account a more recent answer by @MBaz (https://dsp.stackexchange.com/a/60750/41790), there is still one loose end. To see this, consider two Gaussians, $G_1(x)$ and $G_2(x)$, and adding them. Their sum function, i.e., $G_1(x)$ + $G_2(x)$, does not look like either function. Of course, we know the sum function's two constituent summand functions. But the sum function does not show them: it only shows their sum.

Now square the sum function. Then the pointwise squaring works as expected. Alternatively, if we use the expansion, then we get three constituent functions: $G_1^2(x)$, $G_2^2(x)$ and $2G_1(x)G_2(x)$. So the squared sum function is just the sum of these three constituent functions and none of the three are 'visible' to us: we only see the sum of the three constituent functions. So this is why we can never see the 'extra' peak functions. Indeed, we cannot see any of the constituent functions. We can only see the their sum and should not have expected that sum to display its constituent functions.


Why are all real frequencies located on the Imaginary axis of the s-plane. I'd like to understand intuitively why this is true? so my question is the following : "Why are all real frequencies located on the imaginary axis of the s-plane. I'd like to understand intuitively why this is true?"

My guess would be that any "real" frequency has a constant amplitude in the time-domain, which is only true for frequencies with real part equal to zero.

Is this right? Or how should I know intuitively that the statement about real frequencies lying on the imiganiry axis is true?

Take the complex (Laplace transform) variable $s$, and look at the basic exponential signal

$$x(t)=e^{st}\tag{1}$$

With $s=\alpha+j\omega$ we have

$$x(t)=e^{(\alpha+j\omega)t}=e^{\alpha t}e^{j\omega t}=e^{\alpha t}(\cos(\omega t)+j\sin(\omega t))\tag{2}$$

From $(2)$ you see that the real part $\alpha$ takes care of the exponential amplitude, whereas the imaginary part $\omega$ determines the frequency of the signal. If $\alpha=0$ (i.e., $s=j\omega$ is on the imaginary axis) you have an undamped sinusoidal signal with frequency $\omega$.


Applying DTW on MFCC frames I want to use MFCC (Mel-frequency cepstral coefficients) and DTW (dynamic time warping) algorithms to make a speech comparison script. I calculated MFCC  and got 13 coefficients for each frame.

How do I calculate DWT from that? Should I discard the frames and just look at coefficents as an array. I read in this post: Speech recognition using MFCC and DTW(Dynamic Time Warping)? that you should calculate DTW between each frame, but I don't understand what I should do after I have calculated all those DTWs. If I have one audio sequence of $n$ frames and the other one of $m$ frames, then I get $n \times m$ DTW costs. And, where do I go from there.

Well I think I found a solution to this, mostly by trial and error.

So what I did and its giving quite good results is :
1) If we have M(No. of frames) X N(13;No of MFCC Coefficients) matrix;

ie X=X[1 2 .....N(13) ; 1 2.....M(No of frames) ]

2) Find the DTW of [X(1;1 2 3.....M)] with the DATABASE [REF(1;1 2 3...M)]
   Here we are taking K=1(Not 2...12)
   Similarly find DTW of [X(2;1 2 3...M)]  and DATABASE [REF(1;1 2 3...M)]

3) On generalizing, Find DTW of [X(k(kth MFCC coef); 1 2 3...m)];

4)So at the end we will get : K(13) values ;

5) Find the mean of these 13(K) values and let that be your Reference for      deciding;

*Disclaimer: I am not certain if this is what all the papers imply, I just found it by trial and error and it gives accurate results


Possible ways to determine the direction of moving object? I'm trying to determine the moving direction of an object in x-axis only. Figure below shows part of my effort.

However, I got stuck. As you can see from the middle image, the green represent the dot moved to right (positive x), the black represent overlapped area and the purple represent the reference image relative to the moving image.
In bottom image, I got two directions for that one object. Left one is going through positive transition, and right negative transition. I mainly use the intensity differences to determine the direction., i.e. Moving image intensity - reference image intensity.
Is there other possible ways to determine the direction of a moving object (to and fro in x-direction only)? Or is there some ways of correction I can apply to?
Thanks in advance guys!

If your images are as clear (I specify clarity here because it is easy to threshold in this case) as your examples here, one easy way of computing the direction vector of a moving object would be to calculate its center of mass across frames.

For this it's as simple as thresholding first (pixels which are part of the object are equal to 1, all other pixels in the background are 0), then taking the average x-coordinate of these set pixels. Find this average x-value for all your frames, and just take the difference between consecutive frames. You determine whether a positive difference corresponds to movement in the left or right direction.


Ideal sampling - question about the 1/T scaling factor Sources discussing spectrum of sampled signals (under 'hypothetical' IDEAL SAMPLING condition) show that the original message spectrum gets replicated at integer multiples of the sampling frequency.

It's noticed that the sources show that the replica spectra (seen in the sampled spectrum) are scaled by a factor 1/Ts .... or 1/T (where T aka Ts is the sampling period, and 1/Ts is the sampling frequency).

My question I'd like to ask is --- it appears that when the sampling period gets relatively small (or relatively tiny), which means sampling frequency gets relatively large, then the scaling factor 1/T becomes relatively large.

So if the sampling period 'approaches' zero, then the 1/T scaling 'approaches' an infinite value, right? In that case, will each of those peaks of the replica spectra (in the spectrum of the sampled signal) become infinite?

I'm not yet able to understand the implications of the factor 1/T for the case where the sampling period becomes really tiny. It just seems to me that the spectral peaks work towards becoming infinite as T gets mathematically reduced toward zero. Have I understood this correctly? Are the spectral peaks meant to become larger and larger with reduced sampling period (aka larger sampling frequency)? I don't understand why the amplitudes can become really large if T is really small. Will this lead to an infinite power or infinite energy condition?

The main thing I'm currently confused about is - if the peak value of the original signal spectrum is 'A', then the peak value of the sampled signal spectrum is A/T (ie. A/Ts), but since 1/Ts is typically a finite number much greater than 1, then the peak value of A/Ts could be some pretty big number. It's as if there's some kind of huge amplification factor introduced due to 'ideal sampling'. Are we meant to just accept this 'big' 1/T amplification factor?

Thanks all!

Yes you are right.

The impulse modulation description of ideal sampling operation introduces an amplitude scale of $1/T$ on the spectrum of the sampled signal, which is a periodic replica of the original bandlimited signal spectrum.

Note that as the sampling period $T$ gets shorter, the scale factor gets larger and the separation between repeating replicas gets larger too. So in the limit you will be left with a single copy of the original spectrum with infinite amplitude scale.

The mathematical definition of energy in the sampled signal is always infinite, regardless of the sampling period (except possibly at the limit of $T = \infty$).

This is an ideal operation. For the practical applications what you are interested in is the finite valued samples of the signal at the sampling instants; the sample values (envelope) will be the same irrespective of the sampling period.


How to do De-Houghing of a Hough transform'ed Image? I'm working with code found at Rosetta Code for creating a Hough transform. I now want to find all the lines in an image. To do so I need the ρ and θ values of each of the peaks in the Hough space. A sample output for a pentagon looks like this:

How can I find a single [θ,ρ] coordinate for each of the 'hot spots' visible in the Hough space?

You are  finding the coordinates of the peaks and then uses the axis to scale those into  [θ,ρ] coordinates.

Depending on how noisy the data, how many false peaks you expect and how much time you have, there are a few ways of doing it. Easiest is to pick some level that is a a real peak, cut of all data below that and then do a center of gravity on each peak to get it's center.

You could also erode/dialte the image until each peak is a single pixel.


Designing an experiment to show blue noise not adding noise below nyquist and not aliasing above I read this paper from 1983 "Spectral Consequences of Photoreceptor
Sampling in the Rhesus Retina".

This paper talks about blue noise sampling patterns in the retina and says that:

  The results
  (Fig. I) indicate that throughout the retina
  the cones provide a novel form of
  optimal spatial sampling: optimal in the
  sense that minimal noise is introduced
  for spatial frequencies below the nominal
  Nyquist limits implied by local receptor
  densities (the limits that would obtain if
  the cones formed a regular lattice), while
  spatial frequencies above the local Nyquist
  limits are not aliased back into
  conspicuous moire patterns but instead
  are scattered into broadband noise.
  Thus, the visual system avoids the aliasing
  distortion of high frequencies inherent
  in any regular arrangement of image
  sampling elements and simultaneously
  minimizes sampling noise for low frequencies
  that fall within its potential
  Nyquist bandwidths. These advantages
  stem from a quasi-random (that is, random
  but not Poisson) spatial sampling
  scheme that apparently has not been
  used in man-made image-recording devices.

I take that to mean that blue noise:

Doesn't add (much) noise when sampling frequencies below nyquist.
Doesn't alias for frequencies above nyquist.

This is in contrast to:

Uniform sampling which doesn't add noise to any frequency, but aliases frequencies above nyquist.
White noise (uniform random) sampling, which adds noise to all frequencies, but doesn't alias for frequencies above nyquist.

That makes sense, but now I'm trying to think of how I could make a simple experiment to show this being true - preferably in the context of 2d sampling and 2d images (image processing / image sampling), but 1d would be ok too.

I get that I could use blue noise to sample sine waves of varying frequencies but am unsure how exactly that sampling would work.

I know how to generate blue noise sample points in any dimension using Mitchel's best candidate algorithm.

I'm guessing I would generate N sample points and then do a reconstruction filter on those sample points, but thinking about that, I can't see how noise would present itself (although, aliasing would just be "missed frequencies"), and frankly I'm not real sure how to do a decent reconstruction filter.

Am I on the right track for making a decent but simple experiment to show these statements being true empirically?  Or should I be doing something different?

Interesting! Please do share your results with us!

The reconstruction is indeed the difficult bit. Maybe you can look into NFFT (non-uniform FFT) to obtain the Fourier transform of your irregular sampled image and determine its frequency content. You’d expect a clear peak in the case of a sine wave below Nyquist, or an aliased sine wave. Noise would show up as additional (small) local maxima.


How to Generate Band Limited Gaussian White Noise in MATLAB? In some papers, I read that the additive noise is band limited Gaussian white.
How can I simulate this type of noise use MATLAB?

You would generate bandlimited Gaussian noise by first generating white noise, then filtering it to the bandwidth that you desire. As an example:

% design FIR filter to filter noise to half of Nyquist rate
b = fir1(64, 0.5);
% generate Gaussian (normally-distributed) white noise
n = randn(1e4, 1);
% apply to filter to yield bandlimited noise
nb = filter(b,1,n);


Using Parks-Mcclellan (or others) for FIR differentiator? I have never seen the use of Parks-Mcclellan for differentiator - most uses tend to be lowpass filters. Is this because Parks-Mcclellan implementation of FIR differentiator has numerical problems and cannot be successful?

It seems though that Parks-Mcclellan for a differentiator should be more studied and obvious, because in terms of fit of a function, a differentiator has frequency response that is basically a sloped straight line or an identity function.

In general, I am having some difficulty finding established optimal ways to implement FIR differentiators. Are there other known ways that are commonly used for FIR differentiators? What about lowpass differentiators, where frequency response basically is a differentiator up to some frequency and from some frequency to $\pi$ it is filtered out, with transition band in the middle?

There is no fundamental difference between the design of standard frequency selective filters (such as low pass, high pass, etc.) and differentiators. The only difference is that for a (linear phase) differentiator, the desired frequency response must be formulated as a purely imaginary function of frequency, unlike standard frequency selective filters, for which the desired response is purely real-valued. But that's a minor issue, and all important algorithmic details remain the same. Generally, algorithms only become different as soon as the desired response must be defined by a complex-valued function (with non-zero real and imaginary parts).

E.g., the Matlab function firpm, which implements the Parks McClellan method, can design FIR differentiators if you supply differentiator as filter type. The same is true for the least squares filter design function firls.

In sum, designing linear phase FIR differentiators is no more difficult than designing any other FIR filter.

Concerning low pass differentiators, take a look at this answer.


Fourier Transform minimum sample requirement Assuming a signal is sampled adequately, what is the minimum size of an FFT window that would allow detecting a specific frequency?  Is it necessary to have samples for at least one complete period at that frequency?

Please define "specific frequency".

Any signal with no frequency components above the Nyquist limit can be exactly recorded and re-constructed.

However, the number of samples determines the frequency resolution.

if you mean "i want to be able to tell the difference between silence and a 1000hz signal, then you would need 4 samples taken less 1/2000th of a second apart - using that you can solve for the 0 frequency signal (dc offset), and the strength of the other frequency.

if you mean "i want to be able to separate a 900hz signal from a 1000hz", then you need to resolve 11 frequency windows, ie, you'll need 22 samples taken at at least 2000 samples per second; then you'll have measured "frequencies around 0, 100,200,300...900,1000 hz", and you can answer "do i have a 900hz signal".

if you want to measure to the nearest hz the frequency of a N hz signal, then you would need to take at least 2N samples, taken at atleast 2N samples/second.


WOLA vs Simple Hamming Window on signal before FFT Im working on a project in which I should FFT a signal coming into my FPGA from ADC.
I want my FFT to show the peak values as clear as possible without the side lobes.

The solution should be straight forward: Window the time domain signal before FFT and I will get better results in the frequency domain, for example using Hamming Window.

But I encountered another approach to this, apply the WOLA (Weight, OverLap, Add) algorithm on the signal before the FFT.

So I used Matlab to check the results.

Lets assume I am capturing this noisy signal through the FPGA:

so first I will check what is the FFT without any window:

I want to decrease the side lobes of this fft, so the straightforward approach is to create a hamming window, and multiply it with the signal, like this:

% Generate Hamming Window using Equation:
n = [1:1024];
w = 0.54 - 0.46*cos(2*pi*n/1024);
% Apply Hamming Window on Input Signal:
hamming_sim = din_window.*w;

and now lets FFT this new windowed signal:

fft_hamm_sim = abs(fft(hamming_sim));
fft_hamm_sim(1) = 0;
fft_hamm_sim(2) = 0;
fft_hamm_sim(1024) = 0;

OK so I am getting a good result as expected.

I encounterd another approach which requires taking the time domain hamming_window*input_signal, and apply WOLA on it and this is how I did it (I will be glad if someone will alert me if its not the right way):

I took the windowed signal and added the first array element with the middle array element and then the second array element added with the middle+1 array element, like this:

    wola_sim = zeros;
    i = 1;
    for k = 1:512
        wola_sim(i) = hamming_sim(k) + hamming_sim(k+512);
        i = i+1;
    end

and now on this signal i applied FFT:

fft_wola_sim = abs(fft(wola_sim));
fft_wola_sim(1) = 0;

It is still decreasing the side lobes but it is easy to see that the result worse.

My questions:

Am I applying the WOLA algorithm right? Apply hamming window and than add the first half of the array with the second half.
Suppose I apply the WOLA right, what is it used for? Is it really suppose to help me supress the side lobes before the FFT? Why I getting worse result than only hamming window?
For my needs, is it okay to just use the Hamming window?

I'll try to briefly answer your questions.

No, WOLA (weighted overlap and add) is applied directly on the result(s) of an inverse FFT operation. The IFFT results ( time segments ) are first multiplied by a particular window function (e.g Hamming ) and then overlap-added to produce the final time domain signal ( waveform ).
As explained in #1, WOLA is used to generate a time domain signal from multiple, overlapping ifft segments. If you have only a single fft segment (from your FFT analysis), then there's no need for (W)OLA at all.
Yes


blob detection using Laplacian of Gaussian filter I'm learning feature-detectors from this lecture notes, and I don't quite understand the Normalized Laplacian of Gaussian filtered image.

Here is the original image:

this is the output presented in the lecture notes, filtered by Normalized Laplacian of Gaussian with $\sigma=2.502$:

and here is mine, using scipy.ndimage.filters.gaussian_laplace with $\sigma=2.502$:

Well, my output image is quite different from the one in the lecture notes. The background of the lecture's output is totally black, even if I did some thresholding, my output is still different from the lecture's.

Moreover, I noticed that edges almost disappeared in the lecture's output, and around each detected blob, there seems to be a ring surrounding it, like this:

Why is that? Am I using Laplacian of Gaussian wrong?

Most likely your images look different from the ones in the lectures because of scaling. Note that the result of the convolution with a Laplacian filter will have positive and negative values. What the resulting image looks like depends on the data type of the array, and on the range to which the values are scaled.  For example, if you store your filtered image in uint8, without doing any scaling yourself, the negative values will be truncated to 0.

Typically, when an image has negative values, it is scaled so that 0 corresponds to a neutral gray, the minimum value to black, and the maximum to white. It looks like this is what is happening in your code.

On the other hand, in this particular case you are probably interested in both maxima and minima of the Laplacian (dark on bright and bright on dark blobs), so you may as well just take the absolute value of the filtered image.  Then you will not have any negative values, which makes the scaling easier.

As far as your question about the edges, this is exactly what should be happening. The value of the Laplacian at an edge should be close to zero. This property is used for very precise edge localization by finding the zero-crossings of the Laplacian of Gaussian.


(python+opencv) haar cascade training using cloud computing I will train my own  cascade, using python and opencv. I want it to be as precise as possible, so as I read it is common that training takes one or two weeks to compile. I want to use cloud computing to not wait that long. I research but cannot find a single guide that explains what to do. I have pozitive and negative images, and a python code. How I can compile my code on cloud with least additional efford? (till now I just use my desktop, just wrote and compile scripts. I have no experience on cloud computing.) I found zero answers when searching opencv,python,cloud computing together. topics like this are for permanent web applications, I guess. I just want to run the code and get results, basically.

I created a google cloud platform compute engine. It gives you the option to choose the number of CPUs and ram size to choose from. So I experimented with a couple of combinations.

4 CPU 26gb ram - after starting the training I realised it used only ~7 GB of ram with all 4 CPUs running at 100%
8 CPU ~25gb ram - used around 9gb ram
24 CPU ~35gb ram. - used around 9gb ram.

All three instances ran at 100% CPU usage. With the increase in the CPU count, the faster the output files were generated.


Does the OFDM MMSE Channel Estimator Require Prior Channel Knowledge? The OFDM MMSE channel estimator seems to require a-priori knowledge of the channel to calculate the MMSE estimate.  See equation 8 here: http://cache.freescale.com/files/dsp/doc/app_note/AN3059.pdf.  The autocovariance matrix explicitly includes the actual channel H or g.

However, you presumably don't have this information, otherwise why would you be estimating it?
If you are able to deal with the complexity of the MMSE estimator, how do you implement it in practice without prior knowledge of what the channel is?

It's true that the MMSE channel estimator requires some a-priori channel knowledge, namely the autocovariance (or autocorrelation) matrix of the channel and the signal-to-noise ratio (SNR). Note that these are statistical quantities while the aim of channel estimation is to estimate the actual channel transfer function at a certain point in time.

The autocovariance matrix $\underline R_{HH}$ is based on some statistical channel model, e.g. the Rayleigh fading channel or AWGN channel. As you have observed correctly, the channel matrix $H$ is unknown and thus $\underline R_{HH}$ cannot be calculated from $H$. It's rather an assumption based on theroetical models or earlier measurements.

The SNR can usually be estimated from the received signal by measuring the noise power on zero subcarriers or analyzing the received constellation diagram.


Best beginner DSP book for developing plugins I'm a computer science student and I love producing music and I really want to get into developing audio plugins. I'm currently learning the JUCE framework and am fairly comfortable with C++. I'm also good with maths.

I've been doing some research but I'm still not sure as to what book I should use to reference/learn digital signal processing. I've been looking at the DAFX book but its done in MATLAB not C++.

I need a book thats suited for beginners and for developing audio plugins.
Any recommendations?

There are a number of books in addition to Zölzer's DAFX, the one I would recommend is:

The Audio Programming Book by Boulanger...

Also have a look at these:

Designing Sound by Farnell

Designing Audio Effect Plug-Ins in C++ by Pirkle

Designing Software Synthesizer Plug-Ins in C++ by Pirkle


Frequency spectrum of an amplitude modulated signal I am talking about a DSB- Suppressed Carrier amplitude modulation.
Let message signal be $m(t)=A_m \cos{\omega_c t}$ and carrier signal be $c(t)=A_c \cos{\omega_m t}$
Then the amplitude modulated signal is
$$s(t)=m(t)c(t)$$
$$s(t)=A_c A_m \cos{\omega_c t} \cos{\omega_m t}$$
$$s(t)=\frac{A_c A_m}{2} \{\cos{(\omega_c+\omega_m) t} + \cos{(\omega_c-\omega_m) t}\} \tag3$$
In frequency domain $s(t)$ is

From eq 3 it is visible that two frequency components $\omega_c+\omega_m$ , $\omega_c-\omega_m$ are present.
Then in frequency domain $S(\omega)$ why isn't there an impulse at these locations ? (and at their negative frequencies, of course)

With $s$, $m$, and $c$ as you display, then the diagram of $S(\omega)$ is incorrect.  If it were, the spectrum of $m$ would be a triangle centered at zero frequency, which it is not.

$$
S(\omega) = \frac{A_cA_m}{4} \left[\\
\delta\left(\omega - (\omega_c + \omega_m)\right)\\
+ \delta\left(\omega - (\omega_c - \omega_m)\right)\\
+ \delta\left(\omega - (-\omega_c - \omega_m)\right)\\
+ \delta\left(\omega - (-\omega_c + \omega_m)\right)\\
\right]$$


Approximating the frequency response of a median filter A median filter is a non-linear and lossy process, so it doesn't have a closed form frequency response as would a FIR filter (say a box filter of the same length) in an LTI system.

But how closely can something similar to a frequency response of a median filter be approximated?
How would this scale with the length of a median filter?
Under what conditions or for what class of signals might this approximation be ballpark "close"?
For what class of signals might this approximation be very inaccurate?
What kinds of frequency domain distortion or additive noise does a median filter produce?

For a start, any non-linear system will not have an easily-identifiable frequency response. So, it's really a nonsensical question. I intend no offense; nonsensical questions are often the  most enlightening!

However one way to try to answer your question is to assume that the LTI filter involved is the mean (rather than the median) of the windowed data.

Then your question:

  Under what conditions or for what class of signals might this approximation be ballpark "close"?

becomes:

  Under what conditions or for what class of signals might the mean be ballpark "close" to the median.

In that case, for a purely stochastic signal, the mean and median are similar when the probability density function (PDF) of the signal is symmetric about the mean.

  For what class of signals might this approximation be very inaccurate?

When the signal's PDF is "very" asymmetric.


z-space replaces multiplication with convolution and frequency eigenvectors The frequency domain is allegedly preferred because it replaces convolution of complexity $n^2$ with a diagonal matrix multiplication. Yet, I see that in z-domain we have multiplication of polynomials, $(a_0z+a_1z^{-1}+a_2z^{-2}+\ldots)(b_0z+b_1z^{-1}+b_2z^{-2}+\ldots)$, which is, despite its name, is nothing more than convolution with $n^2$ operations: you need multiply one polynomial with every term of the other and combine the terms. This looks like square convolution matrix rather than diagonal one. What is the benefit?

In this regards, I recalled my old question: why nobody speaks about eigenfunctions of the s- and z-domains? Everybody speaks that transfer function is linear and diagonal in the frequency (s- and z-domains). But, when you speak about linear operators, you must bring up their natural basis, especially when speak about diagonalization. How do eigenvectors look like in the Laplace and z-domain? What is the basis of this domain? Might be you do not need to prove the convolution theorem then?

There may be a confusion here about how to cascade filters versus how to actual perform filtering.

There many different ways to implement linear time invariant filtering. If the filter can be represented as an FIR filter than two of the choices are time-domain convolution or frequency domain filtering such as overlap add. The later does indeed just require complex multiplication in the frequency domain (plus the actuall FFTs as well). For large filter length frequency domain processing is a lot more efficient but induces significant latency.

There are many, many other ways to implement filters as well.

Not sure what you mean by "diagonal matrix multiplication". For a system with a single input and a single output, this is just a simple vector multiply and there is no need for a matrix. While you can express time domain convolution as a matrix multiplication, it's not particularly efficient to actually implement it this way, so the matrix representation is mainly used for mathematical analysis.

EDIT: Lets use a specific example

Let's say you have a real signal of 1024 points and a  real FIR filter of 1024 taps.
Doing a direct convolution takes (roughly) 1024*1024 multiply-add operations for a total ca. 2.1 million scalar operations.

Doing it in the frequency domain requires a forward FFT, complex multiplication in the frequency domain and an inverse FFT. A forward FFT of length 1024 takes about 5*1024*log2(1024)= 51200 scalar operations. A complex multiply takes about 6 scalar operations. So the total would be 2*51200+6*1024 or roughly 110k operations. In this particular example the frequency domain multiplication is about 20 times more efficient than the time domain convolution.

DISCLAIMER: This example glosses over a lot of details of the frequency domain method, in particular zero-padding, doing efficiently complex FFTs of real valued signals and overlap handling. However, in terms of computational complexity these mostly cancel each other and don't make a big difference in total operation count.


Filter problems in python I am recording a 30 second 1000Hz Sine wave, saving it to a list, splitting the list into 5 second sublists and I would like to apply a notch filter to each sublist to check when the signal drops out for a few samples (I'm filtering out the 1000Hz frequencies). I am using python along with scipy for signal processing. I had some success using this method when I recorded the signal and applied a filter to the entire signal at once, but whenever I try to apply a filter consecutively to 5 second chunks of the signal, I seem to detect glitches that aren't there.
Here is the code for my filter.

   f0 = 1000.0  # Frequency to be removed from signal (Hz)
    Q = 30.0  # Quality factor
    w0 = f0/(fs/2.0)  # Normalized Frequency
    # Design notch filter
    b, a = signal.iirnotch(w0, Q)
    sos = signal.tf2sos(b,a)
    zi = signal.sosfilt_zi(sos)*lists[0][0]
    for list in lists: #each list is 5 seconds long
        Time=np.linspace(0, len(list)/fs, num=len(list))
        plot_graph(Time,list,fs)
        y_sos, zi = signal.sosfilt(sos, list, zi=zi*list[0])
        plot_graph(Time,y_sos,fs)

The graph below shows the first filtered 5 second chunk.

The large peak just after one second is a genuine dropout but the 3 smaller peaks should not occur at all.

The graph below shows the next filtered 5 second chunk

There are no glitches at all in this section, yet the small peaks still occur. I would also ideally like to get rid of the peak that begins at 0 and 5 seconds of each chunk.
Can anyone tell me what is wrong with my filter/method?

Edit:

Changed my code to the following:

    f0 = 1000.0  # Frequency to be removed from signal (Hz)
    Q = 30.0  # Quality factor
    w0 = f0/(fs/2.0)  # Normalized Frequency
    # Design notch filter
    b, a = signal.iirnotch(w0, Q)
    zi = signal.lfilter_zi(b,a)*lists[0][0]
    for list in lists: #each list is 5 seconds long
        Time=np.linspace(0, len(list)/fs, num=len(list))
        plot_graph(Time,list,fs)
        y, zi = signal.lfilter(b, a, list, zi=zi)
        plot_graph(Time,y,fs)

I'm still obtaining the same problem of strange peaks occuring where they shouldn't although it seems to remove the peak that occurs at the beginning and the end of each chunk. The following graph shows the first filtered 5 seconds, the big peak after 1 second being a genuine dropout, the smaller ones are not sropouts.

The following graph shows the next filtered 5 seconds where there are no dropouts and these peaks should not be occuring at all.

Edit 2:
Zoomed in miniature peak

y_sos, zi = signal.sosfilt(sos, list, zi=zi*list[0])

should just be

y_sos, zi = signal.sosfilt(sos, list, zi=zi)

same as the question you asked yesterday


Successive dynamic time warping with a set of window sizes I am applying DTW to a set of data. I haven't figured out a solution yet, still searching for ideas, and I'm currently trying to apply DTW successively with a monotonic set of windows.

Let's say i have a reference and test signals, which have been both resampled at 10 millis (reference is usually a 10 samples or so of length), the test signal is varying,but 3k of data samples is not unusual. My intention is to run DTW with window size of let's say 10,20,...50 and see how far it gets me. Question is, can I optimize this? Example: starting with 50 and calculating properly the DTW for this, are there shortcuts when progressing down to 40, 30 and so on (or the reverse, starting with 10 and progressing to 50)?

Can the processing effort be reduced under 5 full DTWs by reusing calculations for a different window size, yet applied to the same two signals?

For instance, I reinitialize the cost matrix only once per a set of 10...50 window DTWs. Cost matrix will have a 3K by 3K size after resampling, and reinitialization for constrained DTW, in my understanding, means resetting the elements to an arbitrarily high value. Checking with matlab, yielded the same results as reinitializing every time before DTW, but the running time was reduced to a third. Still, this has been tested insufficiently.

Please note this is not a question about generally optimizing DTW (currently not investigating bounding or anything else), but one with specific application to this use case.

Hmm

There are few properties to help you.

Let Wi be a warping window of width i

It is the case that DTW_wi(A,B) >= DTW_wi+1(A,B)

Other than that, there is little to say.

If you have not already, you should read
http://www.cs.unm.edu/~mueen/DTW.pdf
and
http://www.cs.ucr.edu/~eamonn/CIKM_2016_DTW_clustering.pdf

If you are only interested in the smallest items in the matrix, then using lower bounds (the LB_Keogh) will give you an easy 100-fold speed up.


How to design a digital filter in python that will run over an uC? I am trying to implement a digital filter over a uC (it doesn't really matter which filter and which micro controller because I'm looking forward to learn how to do it in the future with different filters and different microcontrollers).
I've been told that you can design, implement and debug a digital filter in python and when everything is ready you can port the code to C without changing anything. How can I do that? I've been searching for a while and I can not find how to do this.

Also, is there a way to plot the transfer function of an implemented fiter?. I mean the actual filter, a function made by me, that takes the input values and performs the calculations. I don't want to plot the transfer function of a filter made with functions like scipy.signal.butter. I want to plot the transfer function of a filter made with a for, some multiplications and sums.

I would really appreciate any help or information that you can provide me.

Filter representation and design

A DTLTI IIR filter is characterized by its transfer function $ H(z) = \frac{Y(z)}{X(z)} = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2} + \dots + b_{P} z^{-P}}{a_0 + a_1 z^{-1} + a_2 z^{-2} + \dots + a_{Q} z^{-Q}} $.
The transfer function is closely related to the difference equation:
$$ y[n] = \frac{1}{a_0} \left( \sum_{i=0}^{P} b_i x[n-i] - \sum_{j=1}^{Q} a_j y[n-j] \right) $$
($x[n]$ is the input at time step $n$, $y[n]$ is the output at time step $n$).
This is equivalent to
$$ \sum_{i=0}^{P} b_i x[n-i] = \sum_{j=0}^{Q} a_j y[n-j] $$
By applying the Z transform to both sides, using the time shift property, and rearranging some factors, you arrive at the formula for the transfer function $H(z)$ mentioned previously.

When you design a Butterworth filter (or any IIR filter) using SciPy, it returns the $b_i$ and $a_j$ coefficients. These coefficients determine the transfer function and therefore the frequency response of the filter. You can use the freqz function to calculate this response (essentially, it evaluates $H(z)$ along the unit circle $z=e^{j\omega}$).

For example, using Python:

from scipy.signal import butter, freqz
import matplotlib.pyplot as plt
from math import pi
import numpy as np

f_s = 360    # Sample frequency in Hz
f_c = 45     # Cut-off frequency in Hz
order = 4    # Order of the butterworth filter

omega_c = 2 * pi * f_c       # Cut-off angular frequency
omega_c_d = omega_c / f_s    # Normalized cut-off frequency (digital)

# Design the digital Butterworth filter
b, a = butter(order, omega_c_d / pi)
print('Coefficients')
print("b =", b)                           # Print the coefficients
print("a =", a)

w, H = freqz(b, a, 4096)                  # Calculate the frequency response
w *= f_s / (2 * pi)                       # Convert from rad/sample to Hz

# Plot the amplitude response
plt.subplot(2, 1, 1)
plt.suptitle('Bode Plot')
H_dB = 20 * np.log10(abs(H))              # Convert modulus of H to dB
plt.plot(w, H_dB)
plt.ylabel('Magnitude [dB]')
plt.xlim(0, f_s / 2)
plt.ylim(-80, 6)
plt.axvline(f_c, color='red')
plt.axhline(-3, linewidth=0.8, color='black', linestyle=':')

# Plot the phase response
plt.subplot(2, 1, 2)
phi = np.angle(H)                         # Argument of H
phi = np.unwrap(phi)                      # Remove discontinuities
phi *= 180 / pi                           # and convert to degrees
plt.plot(w, phi)
plt.xlabel('Frequency [Hz]')
plt.ylabel('Phase [°]')
plt.xlim(0, f_s / 2)
plt.ylim(-360, 0)
plt.yticks([-360, -270, -180, -90, 0])
plt.axvline(f_c, color='red')

plt.show()

The coefficients can be calculated manually, as explained here (do note that the indices of the coefficients are flipped compared to the formulas above), but it's much easier to use filter design tools like SciPy to calculate them.

Implementation

The difference equation can be used directly to implement the filter. Just loop over the previous inputs and outputs, multiply everything with the respective coefficients, and sum all terms.

This means that once you have designed your filter in Python, you just need to copy the $b$ and $a$ coefficients to your microcontroller to use the filter.

A possible C++ implementation could be:

class IIRFilter {
  public:
    template <size_t B, size_t A>
    IIRFilter(const double (&b)[B], const double (&_a)[A]) : lenB(B), lenA(A-1) {
      x = new double[lenB]();
      y = new double[lenA]();
      coeff_b = new double[2*lenB-1];
      coeff_a = new double[2*lenA-1];
      double a0 = _a[0];
      const double *a = &_a[1];
      for (uint8_t i = 0; i < 2*lenB-1; i++) {
        coeff_b[i] = b[(2*lenB - 1 - i) % lenB] / a0;
      }
      for (uint8_t i = 0; i < 2*lenA-1; i++) {
        coeff_a[i] = a[(2*lenA - 2 - i) % lenA] / a0;
      }
    }
    ~IIRFilter() {
      delete[] x;
      delete[] y;
      delete[] coeff_a;
      delete[] coeff_b;
    }
    double filter(double value) {
      x[i_b] = value;
      double b_terms = 0;
      double *b_shift = &coeff_b[lenB - i_b - 1];
      for (uint8_t i = 0; i < lenB; i++) {
        b_terms += x[i] * b_shift[i];
      }
      double a_terms = 0;
      double *a_shift = &coeff_a[lenA - i_a - 1];
      for (uint8_t i = 0; i < lenA; i++) {
        a_terms += y[i] * a_shift[i];
      }
      double filtered = b_terms - a_terms;
      y[i_a] = filtered;
      i_b++;
      if(i_b == lenB)
        i_b = 0;
      i_a++;
      if(i_a == lenA)
        i_a = 0;
      return filtered;
    }
  private:
    const uint8_t lenB, lenA;
    uint8_t i_b = 0, i_a = 0;
    double *x;
    double *y;
    double *coeff_b;
    double *coeff_a;
};

Circular buffers are used to keep track of the previous inputs and outputs. The coefficients are duplicated in a circular fashion as well, in order to simplify the indices in the loop.
If you're not allowed to use dynamic memory, you could easily create a generic class with stack allocated arrays.

For higher order filters, numerical issues may arise, as a small rounding/quantization error on one of the coefficients affects the location of all of the zeros/poles of the filter.  This is especially problematic for poles that lie very close to the unit circle.
A solution is to factor the transfer function into a product of second order sections. The entire filter is then implemented as a cascade of biquad filters, each implementing one of the sections.
SciPy can produce the coefficients of these sections as well, just pass the output='sos' parameter to the butter function.

Checking the result of the implementation

Numerically calculating the frequency response from a given filter implementation is not straightforward. However, you can compare the impulse response of your implementation with the impulse response produced by SciPy (using scipy.signal.dimpulse).
The transfer function is the Z transform of the impulse response, so if the impulse responses match, the frequency characteristics will match as well.


How to find the time of beginning of the slope of the signal? I have a time signal in MATLAB format. Sampling rate is 100k, so it is roughly 5 seconds of signal. See the link below

https://www.dropbox.com/sh/v7vev2j14i77mlk/AAA-X7JGMAzO6IyElYmFPFwra?dl=0

If you see the freq power vs time graph, you will see something below

Now I want to find the slope of the line on the freq power vs time graph (5k Hz to 15k Hz, 0.5s to 1.5s ) preciously by using the time signal. How should I handle this?

I observe the following about your signal:

of interest are only frequencies below 20 kHz
there's but a single, very dominant tone in that – you've got excellent SNR
The development of frequency over time is either constant, or a linear function of time

So, from that, I'd propose the following steps:

Low pass filter appropriately to restrict your bandwidth to below 15 to 20 kHz
employing a parametric spectrum estimator to find the frequency of your single tone
post-process that frequency estimate to decide whether you're still in the linear or in the constant frequency regime.

Regarding these steps:

might be easiest with just a simple second-order PLL. If you don't know where your tone is going to be at the start of your observation, start with a very high loop filter bandwidth, and as soon as you have a lock, reduce the bandwidth drastically. That second-order control loop actually directly gives you an error signal that you can use for 3.
that's a problem for a bit of low-pass filtering and estimating the frequency/time slope from that, and then checking whether the hypothesis still fits, or the slope has suddenly changed.


Upsampled input to an Adaptive filter? I will try to explain the issue I am having as clearly as possible without going into my coding or maths. I have my own and a MATLAB Central implementation pf standard LMS in MATLAB. Fixed step size. No normalization or other stuff.

I am trying to use it in a system identification setup. I generate a vector of gaussian numbers using "randn" and give the same vector as input and desired response to the LMS filter. Now the estimated weight vector at the end should be a "delta" channel and this is what I get. Then i tried upsampling and interpolating the input vector by an integer number and repeating the same thing. This time around the estimated channel is of the shape of an "Sinc". I gave the interpolated signal as the input and the desired response as before. No changes.

Then i also tried low pass filtering the input vector and repeating the same thing. Again a "Sinc". Has anyone observed this before or know something about this? Please point out my mistake. Any suggestions or a discussion is also welcome.

Prior to upsampling, you have a white signal meaning every single frequency in the Nyquist bandwidth from $-\pi$ to $\pi$ is represented. This is a requirement to obtain an impulse (because the Fourier transform of an impulse is a white spectrum). The reason that white signals are often used as inputs for purposes of system identification is that they excite all frequencies of the system. This is important when trying to determine the system's response.

After upsampling and interpolating by $N$, you obtain a sequence with a spectrum that is more or less white in the bandwidth from $-\frac{\pi}{2N}$ to $\frac{\pi}{2N}$. The rest of the band (if your interpolation filter is working correctly) will be nearly zero or null. That means that these frequencies are no longer contributing to your coefficient calculation. In turn, this means that you cannot obtain a delta function since it requires every frequency present to be represented properly. The closest approximation to a delta function that you can obtain is a sinc pulse. This is why you are observing this result.

Summary: You cannot estimate the response at frequencies that are not contained in your input signal. Your input signal and desired signal are no longer white (because they are bandlimited by the interpolation filter). Therefore, the adaptive filter does not estimate those frequencies. This leads to the sinc response you are observing which is the closest band-limited response to the true system response (a delta function).


Does it make sense to use a wavelet that is equal to a sine of one period? Would I get same results with STFT method (short-time Fourier transform) if I try this?

With all due caution, no in both cases (title and body question). I'll start with the second one.

Continuous wavelets use all dilations of the mother wavelet, which are not accessed with the STFT
The STFT is complex in general, and the windowed sine is not.

For the first one: I  never tried it, and do not remember having seen it in use, and one should check first whether it satisfies the admissibility condition (first thought: seems OK on 0 in Fourier, and decays at infinity, I'd expect families of synthesis wavelets). You can compute the effects in Fourier using SE.Math: Fourier transform of 1 cycle of sine wave.

However, the sine falls abruptly at its edges, without differentiability, in other words it is not smooth enough. So I would not expect a  more interesting interpretation  of the scalogram than with other continuous wavelet.

But it could be used as an exploratory tool (analysis) for some specific signals. You can think about  hybrid methods as well: try to fit (in more or less robust ways) a parametric sine model, at different scales.


Averaging Different Length Speech Signal I have around 10 wav audio files all containing same word spoken i.e. 'Hello' but all having different length.
So lets say 1st file when loaded consist of 20000 samples, 2nd of 15000, 3rd of 18000 and so on they vary. Now I want to create a average template made of averaging all these 10 files, so that I can have better speech recognition result. But the problem is all the files have different length.

How can you average 10 different length signals?

Dynamic Time Warping is pretty well explained on this site. I'll use some of the diagrams from the PPT on that site to explain.

The idea is to divide the signals into segments (frames) and then compare frames sequentially through each signal.  As illustrated below, motion from a segment in one signal to the next segment depends on the similarity to the current segment in the other signal.

For your example of audio files, you usually have features that you extract from them (LPCs or MFCCs for example).  Extract these from each frame.  Use these features as the comparison features used in the DTW algorithm.

  1) Suppose I have 10-10 frames of my 2 signals and each frame has 13 MFCC coefficients so the 2 DTW signals will be all coefficients arranged sequentially?

Well, before DTW you'll have $MFCC_s^f$ for signal $s$ for each frame, $f$:

Signal 1: $MFCC_1^1$,$MFCC_1^2$, $\ldots MFCC_1^{N_1}$

Signal 2: $MFCC_2^1,MFCC_2^2, \ldots MFCC_2^{N_2} $

where $N_1 \not=N_2$ and $N_s$ is the number of frames for each signal, $s = 1,2$.

After DTW, you'll have the best matching path over $N_\mbox{min} = \min(N_1,N_2)$ frames:

Signal 1: $MFCC_1^1$,$MFCC_1^1$, $\ldots MFCC_1^{N_?}$

Signal 2: $MFCC_2^1,MFCC_2^2, \ldots MFCC_2^{N_?} $

Then the comparison will be done by running all $those$ MFCCs sequentially.

  2) How this shortest path can be used to decrease or increase the size of one signal so that they are of same size?

See above explanation: After DTW you'll only have the same number of MFCCs for each and every signal.  All other (unused) frames with MFCCs are discarded for comparison purposes.


How to interpret output of matched filter with complex input? I have implemented a matched filter based on the Fourier Transform approach. In the real numbers domain that means that I use as the coefficients of my filter (B) the inverted time-samples of the signal that I'm trying to find and I compute:

real_output = IFFT(FFT(A).*FFT(B))

Where, if we assume the moment of exact match, A is an input signal of real value samples and B is the signal made from inverted time-samples of A (or I can also conjugate the spectrum of B and I get the required time inversion because the signal is real). In order to avoid aliasing I perform both FFT with length(A)+length(B)-1 samples.
The output of this filter is always real and its maximum, points at the moment in which a match has occurred.

Now assume that instead of a real valued A we have a complex valued A. My questions are:

Is it possible to apply a matched filter for complex domain signals?
Would I have to invert the complex valued A to get the coefficients of my filter or should I do something else like conjugating them?
How do I interpret the output of the filter, should I say a match happens if the absolute value of the output is maximum or just when the real part is maximized?
In some experiments I have noticed that at the moment of a match the output of the filter is real if spectrum(B) is the complex conjugate of spectrum(A), is this correct?

Yes, it is possible (at least on paper or code, since complex signal don't exist physically) to apply a matched filter to complex signals.

This is one way to look at it that I think is illustrative. Assume that a pulse $p(t)$ is real and has energy 1, that is, $$\int_{-\infty}^\infty p^2(t) \, dt=1.$$ The filter matched to $p(t)$ has impulse response $p(-t)$. Then, if the filter's input is $ap(t)$ for real $a$, the matched filter's output (sampled at the appropriate time) is equal to $$ap(t)\star p(-t)=a.$$

If the pulse amplitude is complex, say $a+jb$ for real $a$ and $b$, then we can use the linearity of the convolution operation to calculate
\begin{align}
(a+jb)p(t) \star p(-t) &= \left(ap(t)+jbp(t)\right) \star p(-t) \\
&=ap(t)\star p(-t) + jbp(t)\star p(-t) \\
&= a+jb.
\end{align}

One way to interpret the result is that the filter is matching the real part and the imaginary part independently (but "at the same time"). Recall that $j=\sqrt{-1}$ is orthogonal to the real line; this means that $a$ and $b$ don't "interfere" with each other during the convolution.

To summarize, the answers to your questions are:

Yes.
If you can define the system as using a real pulse $p(t)$, then the matched filter is $p(-t)$. If you must use a complex pulse, then the matched filter is $p^\star(-t)$.
A "match" can be declared when the filter's ouput is close enough (in Euclidean distance sense) to the matched signal (see https://en.wikipedia.org/wiki/Maximum_a_posteriori_estimation for more details).
The output of the filter should be real only when $b$ is zero.


why does my integrator perform noise shaping? I am designing (in matlab) a simple system with a DC signal at 0 Hz and some added noise.
I am performing a differentiation of the signal (by simply subtracting the last values from each other), next i am applying rolling average algorithm, to filter the high frequencies, and then i simply integrate over a number of samples. Using the pwelch() function on the output i noticed that there is some noise shaping going on, and i don't really now where it is coming from...

I am aware of delta-sigma modulators, but this is not a thing i am doing here: for my integrator i am using a very simple algorithm that adds $n$ last samples together. I would expect a simple integrator behaviour here (with frequency response falling 20db/dec).

Why does the noise shaping take place?
https://i.stack.imgur.com/dTfLl.png

I am posting the code for the signal generation, differentiator and integrator below.

clear
home
close all

%generate the signal
    fs = 1;
    T = 1/fs;
    t = 0:T:(2^14);
    pos = 3 * t;
    noise = 2 * randn(1, length(t));
    pos = pos+noise;

    %differentiate the signal
    diff_time = 1;
    diff_out = [];
    len = length(pos);
    for j = 1+diff_time : len
        diff_out(j-diff_time) = pos(j)-pos(j-diff_time);
    end

    %filter the hf with MA filter with 16 taps
    MA_avg = zeros(1,16);
    for i = 1 : length(diff_out)-16
       avg_val = MA_avg(end);
       old_val = diff_out(i);
       current_val = diff_out(i+16);
       MA_avg(end+1) = current_val-old_val+avg_val;
    end

    %integrate
    pos_intgr = integrate(MA_avg, 16);

    %get the spectral power
    NFFT = length(pos_intgr);
    [P, F] = pwelch(pos_intgr,ones(NFFT,1),0,NFFT,fs,'power');
    PdBW = 10*log10(P);
    plot(F,PdBW)
    title("pwelch")
    xlabel('Frequency')
    ylabel('Power spectrum (dbW)')

    function [acu_down] = integrate(signal, taps)
       len = length(signal);
       acu_down = [0];
       i = 1;
       k = 1;
       while i <= len-taps
          sum = 0; %reset the sum
          for j = 0 : taps-1
               sum = sum + signal(i+j);
          end
          i = i+taps;
          acu_down(k) = sum;
          k = k+1;
      end
      %cut out the LSB
      acu_down = acu_down ./ taps;
   end

The dominant effect you see is from the differentiation which is a high pass function.

Observe the spectrum after differentiation alone:

As would be clearer on a log log plot the signal is going up at rate f versus frequency. (consistent with a differentiation).

Next observe the spectrum at the moving average output:

This is consistent with expected: The moving average approximates a Sinc filter response, and the nulls of a Sinc function go down at rate f versus frequency. With the differentiator going up at rate f and the moving average going down at rate f we see the flat response as observed with the nulls at $1/T$ where $T$ is the duration of the moving average.

Your final operation that you call an integrator is actually another average combined with a down-sample (decimation), and with that would result in the low frequency portion of your signal after another 16 point average; so another Sinc filter response of similar form and then the decimated (close in) spectrum out to half of the first main lobe:


Comparing spectrogram similarity over time I have a data set of spectrograms from moving audio recorders. As the recorders get closer, I expect the spectrograms to become more "similar". I expect this "similarity" to be a function of distance and frequency.
I would like to create a 3d plot or image that has distance and frequency on the x and y, and "similarity" on z. I am not sure how to define this similarity. Is there a standardized way of comparing similar signals in this way?

For spectrogram similarity, I suggest using a metric from image processing known as SSIM - Structural Similarity Index Metric. It is a quantity between $[0,1]$ which denotes the similarity between two images.
The paper is called Image Quality Assessment: From Error Visibility to Structural Similarity. The MATLAB function ssimval() will give you the metric.


Compare image segmentation approaches I've got a question related to a comparison of image segmentation algorithms.
I'd like to compare two algorithm results to each other in step one.
In another step I'd like to compare the two algorithms results against a few ground truth segmented images.
Unfortunately I got only a few ground truth data results in form of single images.
So I can't compare all possible images and parts of an image to the two algorithms.

My first idea on step #1 was that I go down on subpixel area and look up which algorithm is more precise in relation to the original image and discuss it without the ground truth data.
I'll take the difference of the pixels on which the difference lasts. If it's too obvious, e.g. if one algorithm fails completely than the other wins automatically.

In the second approach where I've got ground truth data I would measure the absolute difference in pixels of each algorithm on significant parts of an image.

So, what I wanted to know is:
1) Are this correct approaches for comparing segmentation algorithms?
2) How do I measure the difference in step 1? Do I have to measure all segments of each algorithm?
3) What is a good quality criteria? Is it possible to get a quantitative criteria for the two algorithms?
4) Do I have to compare the whole segmented parts of each algorithm to the ground truth data to gain absolutely correct results?

Please let me know if something is unclear or not correctly formulated.
Any help, advice or feedback is appreciated!
Thank you guys in advance.

Sincerely

Typically to compare segmentation algorithms, people use an NPR index. NPR stands for Normalized Probabilistic Rand.

An over segmented or an under-segmented image will typically have negative values of NPR whereas a properly segmented image will have positive values of NPR index.

This type of comparison is typically carried out for segmentation algorithms of the split-and-merge type.

Following this paper and the further publications of the group would be advantageous.


Impulse response of a certain system I know that generally the impulse response of a system is the output when input is an impulse. It is also the differentiation of a step response. So if we have a system described by thi differential equation:
$$\frac{d^2 y(t)}{dt^2}+5\frac{\text{d}y(t)}{\text{d}t}
+6y(t)=-\frac{\text{d}x(t)}{\text{d}t}$$
where $x(t)$ is the input, $y(t)$ is the output, how do I find the impulse response?

My first idea was to let $x(t)$ be a ramp function so that it's derivative would be a step function then I could easily solve the differential equation and differentiate what I got to get the step response then differentiate it again to get the impulse response. But now I'm doubting this path as the ramp function is not a bounded function so the response might be unstable. I'm not sure how to approach this.

Let me provide a method, applied only for finding the impulse response $h(t)$ of an LTI system characterised by an LCCDE of the form $ \sum_{k=0}^{N}{ a_k  {{d^k y(t)}\over {dt^k}}} = \sum_{k=0}^{M}{ b_k  {{d^k x(t)}\over {dt^k}}}$ by using the classical time domain approach. A method which is ignored as it's replaced by transform domain techniques instead...

Note that the input $x(t) = \delta(t)$ is formally a problematic signal (function) and therefore mathematically oriented classical books on ODEs, do generally avoid any discussion of such functions and their solutions, unless the scope of the book specifically includes generalised functions (distributions), usually to be used in some engineering or physical fields.

The procedure: Consider an LTI system which is causal with initial rest conditions. We say that when the input to this system is $x(t)=\delta(t)$ an impulse, then its output $y(t)=h(t)$ is the impulse-response of the system. We find this solution $h(t)$ in two steps, by breaking the LCCDE into two parts as inspired by a serial (cascade) implementation of two LTI systems corresponding to the right and left sides of the LCCDE.

I/O relationship of Part-I is given by the LCCDE: $$\sum_{k=0}^{N}{ a_k  {{d^k y(t)}\over {dt^k}}} = x(t)$$ which represents the first stage with its input $x(t)=\delta(t)$ and we denote its solution (stage1-output) as $h_0(t)$ which is actually the impulse response of Part-I.

And the I/O relationship of the Part-II is given by the equation: $$ y(t) = \sum_{k=0}^{M}{b_k{{d^k x(t)}\over {dt^k}}} $$ which requires nothing but simple summation of its input, $x(t) = h_0(t)$ ,and its derivatives to compute the output as $$h(t) = \sum_{k=0}^{M}{b_k {{d^k h_0(t)}\over {dt^k}}}$$ Therefore we need to find $h_0(t)$ of the Part-I to simply compute the impulse response $h(t)$ of the complete system.

To find the solution of part-I: $\sum_{k=0}^{N}{ a_k  {{d^k y(t)}\over {dt^k}}} = x(t)$ whose solution $y(t)$ is what we call $h_0(t)$ , we we make the following observation that $y(t)$ is composed of two parts as $$y(t) = y_{h}(t) + y_{p}(t)$$ where $y_{h}(t)$ is the homogeneous solution corresponding to $x(t)=0$ for all t, and $y_{p}(t)$ is the particular solution corresponding to input $x(t)$ for $ t > 0$, in particular as what remains as input when time goes to infinity.

Now, we can see that for the input $x(t)=\delta(t)$ the particular solution $y_p(t)$ for $t>0$ is zero as the input is zero for $t > 0$. Then the output simplifies to $h_0(t) = y(t) = y_h(t)$

Now, we have to find the homogeneous part of $y(t)$ as the solution of the equation $$\sum_{k=0}^{N}{ a_k {{d^k y(t)}\over {dt^k}}} = 0$$
following the usual procedure by finding the roots of its characteristic polynomial in complex-s: $$p(s) = \sum_{k=0}^{N}{a_k s^k}$$

which is obtained when we simply insert $y_h(t) = Ke^{st}$ as a solution canditade to the homogeneous differential equation $\sum_{k=0}^{N}{ a_k {{d^k y(t)}\over {dt^k}}} = 0$. The solution is then hypothesised to be $$y_h(t) = \sum_{k=0}^{N}{ A_k e^{st}}$$ for the case when all the roots are distinct. If they are not, you know what to do right ;)
The solution therefore now requires to determine the N unknown coefficients $A_k$ from N-inital conditions, at time t=0 for the framework of this method, possesed by the output $y(t)$ and its derivates up to order N-1.

Remember that this system was defined to be under inital rest with complete zero inital conditions prior to the application of any excitation. Hence if we were to solve the homogeneous solution $y_h(t)$ for an ordinary (formally valid) input signal $x(t)$ to the system, the homogeneous solution will be identically zero, as the initial conditions would all be. (all the coefficients $A_k$ of the exponential terms turn out to be zero). But the magic of the generalised input $\delta(t)$ is that it will set one of those inital conditions to non-zero, thereby, enabling a non-zero homogeneous response to exist, which will become the solution of the part-I as well.

It can be shown that the inital conditions of the part-I due to the excitation $\delta(t)$ will be set as $y(t)=0 ~,~ y(t)' = 0 ~,~ y(t)''=0 ~,~... ~,~ y(t)^{N-1}=1/a_N$ all at $t=0$.

From this information you can proceed to solve the algebraic equaitons as usual and find $h_0(t)$ from which you can compute the impulse response $h(t)$ for all $t > 0$. What happens for all $t <0 $? As we stated, this system was LTI and causal with initial rest condition. For such a system if the input $x(t)$ is zero for all $t<0$ then that output is also zero for all $t<0$ which is achieved by multiplying $h_0(t)$ via $u(t)$ to represent the solution of Part-I for all t.

It is instructive to apply this procedure to your example:
$$y'' + 5 y' + 6y = -x'$$

Lets first solve Part-I:  $$y'' + 5 y' + 6y = x$$
To find $h_0(t)$ when the input is $x(t)$ = $\delta(t)$ first solve the characteristic equation: $s^2 + 5s + 6 = 0$ which yields $s=-2$ and $s=-3$.

The proposed homogeneous solution is therefore $$y_h(t) = A_1 e^{-2t} + A_2 e^{-3t}$$

For which we need 2 inital conditions to find $A_1$ and $A_2$.
These initial conditions are set by the impulse $\delta(t)$ to be: $y(0)=0$ and $y(0)'=1/1=1$. Using these initial conditions we find $A_1$ and $A_2$ to be $A_1=1$ and $A_2=-1$ and this yields the solution to part-I as
$$ h_0(t) = [ e^{-2t} - e^{-3t}] $$ for all $t > 0$.

Also note that this system was LTI and causal with initial rest and input $x(t)=\delta(t)$ is zero for all $t<0$ implies that output is also zero for all $t<0$ which is simplified by using a unit step $u(t)$ function to format the output as:
$$ h_0(t) = [ e^{-2t} - e^{-3t}]u(t)$$ for all t.

Finally, to solve part-II and to get the impulse response of the overall system, apply this $h_0(t)$ as an input to Part-II whose output is simply

$$y(t)= -x(t)' \longrightarrow h(t)= -h_0(t)'$$

$$h(t) = -[ -2 e^{-2t} + 3 e^{-3t}]u(t) - [e^{-2t} - e^{-3t}]\delta(t) $$
$$h(t) = -[ -2 e^{-2t} + 3 e^{-3t}]u(t) - [1 - 1]\delta(t) $$
$$h(t) = -[ -2 e^{-2t} + 3 e^{-3t}]u(t) - 0\delta(t) $$
$$h(t) = [2 e^{-2t} - 3 e^{-3t}]u(t)$$

Note that the impulse at the output is cancelled by the property $f(t) \delta(t-a) = f(a)\delta(t-a)$


Real-time DAQ (Data AcQuisition) device and filtering of voltage signals in audio frequencies I'm working on a frequency mixer project in which two analog voltage signals at two different frequencies (each <5 kHz) are considered. I'd like to add the two fundamental frequencies and return an analog output in real time. Can I do this digitally and how difficult it is compared to an analog approach?

My idea so far includes an ADC with a fast enough sampling rate, Python scipy (for filtering), and DAC for output.

I'm trying to go low cost (<100 bucks). As far as I'm concerned, I don't see a potential problem in this idea as an ADC with Msps can be obtained, the math process is fast enough, and the DAC can also have Msps sampling rate.

But as a chemistry student, I'm no professional at DSP, so please point out some of the potential flaws/drawbacks of this method and tell me if you'd recommend doing this in analog?

Thanks in advance.

To answer my own question here for anyone who would like to do something similar, I ended up designing an analog circuit instead of DSP. I used a signal multiplier chip called AD633 to multiply the two frequencies. Then I use a bandpass filter to extract the sum frequency.


How to convert a sound signal function to wave file? I am new to understanding sound signals. A sound signal, as far as I know, is a real function of time $x:[0,\infty)\to \Bbb R$. For example, $x(t)=\sin(t)$ is a sinusoidal signal.

What is $x(t)$. For example, what are the numbers $x(1)$ or $x(2)$? What do these numbers represent?
How can I produce a sound file e.g. sound.wav from a signal function, e.g. $x(t)=\sin(t)$. Can you give me a simple code (in c#, c++,..) for converting $x(t)$ to a sound file? (I need the code just to see and understand the details).

i don't consider this a "bad" question.  But there is a lot that nano needs to deal with.

first, you must be able to think about conceptually and mathematically converting your continuous-time signal

$$ x(t) = \sin(\Omega t) $$

into a discrete-time signal

$$ x[n] = \sin(\omega n) $$

how $n$ is related to $t$ and how $\omega$ is related to $\Omega$.  this requires knowledge of the sample rate.

finally, after creating an array of samples, you can write that array to a .wav file using a library such as this one from Erik de Castro Lopo.  you might have to massage the code to make it work for C#.

or you can learn about the format of the particular .wav file (nowadays they can be float or int format) and write that file using the C# counterparts to fopen(), fwrite(), and fclose().  but i dunno shit about C#.


Interpolating the peak of a cross-correlation using a centroid Forgive me if this question is poorly worded, I'm not sure if centroid is the word to use here.

Say I want to interpolate the peak of a cross-correlation function in order to get sub-sample delays. I can fit a parabola with three points, or simply upsample my signal through zero-padding/sinc-interpolation. However, I haven't seen any comments (or found any literature) on using some sort of weighted average on the main peak.

In case it's unclear, i mean performing the following:

Finding the maximum of the cross-correlation in order to get a coarse estimate.
$$
\hat{D} = argmax(R_{12})
$$
Taking the P nearest neighbors (where [-P,P] corresponds to the central peak area) and doing the following fine estimate

$$
\hat{d} = \frac{\sum_{k =-P}^P kR_{12}(k)}{\sum_{k=-P}^PR_{12}(k)}
$$

I've tested it for a gaussian distributed, infinite snr signal, critically sampled (so peak is only three samples. Centroid, (yellow) versus some sinc interpolation (purple) and parabolic fitting (orange), and it seems to offer good performance.

It appears to have less maximum bias, despite the discontinuity.
I also tested it with noise and the results were also okay. I get that this requires knowledge of the signal A.C.F. beforehand whereas fitting a parabola doesn't, but is there any other reason why this method is not seen in the literature?

Either way, I'm finding it somewhat advantageous in my application (also not sure wh), so I would like to read some more formal studies or descriptions of the technique. My background is not signal processing, so I'm afraid I'm overlooking some severe shortcoming.

NOTE: I also tested the mean-square error. It seems to perform much worse there,
 but as I said, for my use case it seems to perform well. . (The black line is the cramer-rao lower bound. Not sure why the parabola goes below it, probably because it's a biased estimator...)

EDIT: Signal example of the application (this picture has quite a lot of distortion of the delayed signal, but it's what I have at hand currently).

There are a lot of papers on bin interpolation and some of them have obvious connections to interpolating a parabola and others aren't' so obvious. Given what you presented, I'm sure that you can make the connection between cross correlation and spectral analysis. I suggest  you look at:

M. D. Macleod, "Fast nearly ML estimation of the parameters of real or complex single tones or resolved multiple tones," in IEEE Transactions on Signal Processing, vol. 46, no. 1, pp. 141-148, Jan 1998.
doi: 10.1109/78.651200
Abstract: This paper presents new computationally efficient algorithms for estimating the parameters (frequency, amplitude, and phase) of one or more real tones (sinusoids) or complex tones (cisoids) in noise from a block of N uniformly spaced samples. The first algorithm is an interpolator that uses the peak sample in the discrete Fourier spectrum (DFS) of the data and its two neighbors. We derive Cramer-Rao bounds (CRBs) for such interpolators and show that they are very close to the CRB's for the maximum likelihood (ML) estimator. The new algorithm almost reaches these bounds. A second algorithm uses the five DFS samples centered on the peak to produce estimates even closer to ML. Enhancements are presented that maintain nearly ML performance for small values of N. For multiple complex tones with frequency separations of at least 4π/N rad/sample, unbiased estimates are obtained by incorporating the new single-tone estimators into an iterative “cyclic descent” algorithm, which is a computationally cheap nonlinear optimization. Single or multiple real tones are handled in the same way. The new algorithms are immune to nonzero mean signals and (provided N is large) remain near-optimal in colored and non-Gaussian noise
keywords: {Fourier analysis;Gaussian noise;amplitude estimation;frequency estimation;harmonic analysis;interpolation;iterative methods;maximum likelihood estimation;optimisation;phase estimation;signal sampling;spectral analysis;white noise;AWGN;Cramer-Rao bounds;FFT;amplitude;cisoids;colored noise;complex single tones;discrete Fourier spectrum;fast nearly ML estimation;frequency;frequency separations;interpolator;iterative cyclic descent algorithm;noise;nonGaussian noise;nonlinear optimization;nonzero mean signals;parameter estimation;peak sample;phase;real single tones;resolved multiple tones;sinusoids;unbiased estimates;uniformly spaced samples;Amplitude estimation;Computational efficiency;Frequency estimation;Iterative algorithms;Maximum likelihood estimation;Noise level;Parameter estimation;Phase estimation;Phase noise;Signal processing algorithms},
URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=651200&isnumber=14197

and if you can get into IEEE Explore, look at the papers that cite it. These are the State-OF-The-Art bin interpolation schemes, which you should be making your comparisons to.

The only comment that I have about you centroid technique, would be that the accuracy would seem to depend on there being only a single delay in your data.  My experience is that where there is one delay, there are often many.  Having too wide a window would seem to be more susceptible to additional delays.


What is the advantages of BPSK+OFDM over usual BPSK (with no OFDM Modulation)? The other days I read a useful article (on BER v.s. Eb/N0 curve of ofdm) whose link is as bellow. While watching the plot inside, a question was arisen in my mind. I wonder if the BER v.s. Eb/N0 curve of BPSK alone (without using OFDM) and one using OFDM could be the same in both AWGN and fading (flat and selective) channels? If so, what’s the use of applying OFDM to BPSK alone (expect the simpler equalizer for OFDM case)?

http://www.dsplog.com/2008/08/26/ofdm-rayleigh-channel-ber-bpsk/

Best Regards

In a pure AWGN channel there is no benefit from OFDM because the noise will affect the received signal in the same way regardless of whether sent as a higher bandwidth serial stream or lower rate parallel streams. Similarly, there isn't much benefit in a flat fading environment either because, by definition, each constituent stream of the OFDM signal will be affected by the same amplitude change and linear phase shift.

The main benefit from OFDM comes with frequency selective fading. That is because the effects of frequency selective fading worsen as signal bandwidth increases. With OFDM, each parallel low rate signal experiences a flat fading channel instead a frequency selective channel at the full bandwidth - i.e. each individual channel can be represented by a single complex value instead a multi-tap filter. That means frequency selective equalization is easily addressed, which allows us to scale to much higher bandwidths - and everybody wants higher bandwidth.

Note that the question refers to BER only, so other OFDM related characteristics such as spectrum mask or resource allocation are left out.


How to prove a train of sinc pulses in digital communicaton system are orthogonal to each other? Consider a train of sinc pulses:
$$\phi_n(t)= \frac{\sin(\omega_M(t-nT_s))}{\omega_M(t-nT_s)}\quad; n=0,\pm1,\pm2,\dots$$
$\quad$where,$\quad T_s=\frac{\pi}{\omega_M}$
Now ,in order to show sinc pulses are orthogonal we need to prove:
$$\int_{-\infty}^{\infty}\phi_n(t)\phi_k(t)dt=T_s \delta_{nk} \quad \dots(1)$$
where, $\delta_{nk}$ is kronecker's delta.
So, i began doing it as follows:
$$T_s=\frac{\pi}{\omega_M}=\frac{\pi}{2\pi f_M}=\frac{1}{2f_M}=\frac{1}{f_N} \quad \dots(2)$$
where, $f_N$ is the Nyquist frequency
$$\phi_0(t)=\frac{\sin(\omega_Mt)}{\omega_Mt}=\frac{\sin(2\pi f_Mt)}{2\pi f_Mt}=sinc(2f_Mt)=sinc(f_Nt) \quad \dots(3)$$
Now, $\mathscr{F}\{ sinc(f_Nt) \}=\frac{1}{f_N} rect(\frac{f}{f_N})$ ,where $rect$ is a rectangular function centred at origin and having width= $f_N$
$$\implies \mathscr{F}\{ sinc(f_N(t-nT_s)) \}=\frac{1}{f_N} \exp(-i2\pi f n T_s) rect(\frac{f}{f_N}) \quad \dots(4)$$
Now we can write:
$$\int_{-\infty}^{\infty}\phi_n(t)\phi_k(t)dt=\int_{-\infty}^{\infty} \{ \Phi_n(f) \circledast \Phi_k(f) \} df$$
$$=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty}\Phi_n(\tau) \Phi_k(f-\tau) d\tau df$$
$$=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \exp(-i2\pi \tau n T_s)
\frac{1}{f_N} rect(\frac{\tau}{f_N}) \exp(-i2\pi (f-\tau) k T_s) \frac{1}{f_N} rect(\frac{f-\tau}{f_N}) d\tau df$$
$$=\int_{-\infty}^{\infty} \exp(-i2\pi \tau n T_s) \frac{1}{f_N} rect(\frac{\tau}{f_N}) \exp(i2\pi \tau k T_s) \{ \int_{-\infty}^{\infty} \exp(-i2\pi f k T_s) \frac{1}{f_N} rect(\frac{f-\tau}{f_N}) df \} d\tau \quad \dots(5)$$
The inner integral of $(5)$ can be simplified as:
$$\int_{\tau -\frac{f_N}{2}}^{\tau +\frac{f_N}{2}} \frac{1}{f_N} \exp(-i2\pi f k T_s) df$$
$$=\frac{\exp(-i2\pi \tau k T_s) \sin(\pi k)}{\pi k} \quad \dots(6)$$
So, $(5)$ can be rewritten as:
$$\int_{-\infty}^{\infty} \exp(-i2\pi \tau n T_s) \frac{1}{f_N} rect(\frac{\tau}{f_N}) \exp(i2\pi \tau k T_s) \frac{\exp(-i2\pi \tau k T_s) \sin(\pi k)}{\pi k} d\tau$$
$$=\frac{\sin(\pi k)}{\pi k} \int_{-\infty}^{\infty} \exp(-i2\pi \tau n T_s) \frac{1}{f_N} rect(\frac{\tau}{f_N})d\tau $$
$$=\frac{\sin(\pi k)}{\pi k} \frac{\sin(\pi n)}{\pi n} \quad \dots(7)$$
Now, $(7)$ is even equal to $0$ when $k=2$ and $n=2$
So,where i missed ? any help or suggestions please...

The idea of solving the integral in the frequency domain is good, but you made a mistake rewriting the integral. Note that
$$\int_{-\infty}^{\infty}\phi_n(t)\phi_k(t)dt\tag{1}$$
equals the Fourier transform of $\phi_n(t)\phi_k(t)$ evaluated at $f=0$. As you know, that Fourier transform is given by the convolution of the two individual Fourier transforms of $\phi_n(t)$ and $\phi_k(t)$, respectively:
$$\mathcal{F}\big\{\phi_n(t)\phi_k(t)\big\}=\int_{-\infty}^{\infty}\Phi_n(\xi)\Phi_k(f-\xi)d\xi\tag{2}$$
Evaluating $(2)$ at $f=0$ gives
$$\int_{-\infty}^{\infty}\phi_n(t)\phi_k(t)dt=\int_{-\infty}^{\infty}\Phi_n(\xi)\Phi_k(-\xi)d\xi=\int_{-\infty}^{\infty}\Phi_n(\xi)\Phi_k^*(\xi)d\xi\tag{3}$$
where the last equality in $(3)$ is true because $\phi_k(t)$ is real-valued. Eq. $(3)$ is just Parseval's theorem.
I'm sure you can continue from here and show that the right-hand side of $(3)$ equals zero for $n\neq k$.

Note that the integral you tried to compute equals the inverse Fourier transform of $(\Phi_n\star\Phi_k)(f)$ evaluated at $t=0$, i.e., it equals $\phi_n(0)\phi_k(0)$ which also satisfies
$$\phi_n(0)\phi_k(0)=\delta[n-k]\tag{4}$$


Envelope analyis basic question I have a vibration signal sampled at 25.6 kHz, and I want to perform Envelope Analysis on it.

"Vibration-Based Condition Monitoring of Wind Turbines" by T. Barszcz has the following figure:

I understand why the frequency is changing in the FFT plot during the HP and LP filters, but why is the frequency range changing from (b) -> (c) [0 -> 8 kHz] when performing the abs (rectifier) step?
Cheers!

Edit:
Here is an FFT for "Gearbox High Speed shaft" from my data.

A high-pass filtered signal has lost its low frequencies, so you can refocus the spectrum on a higher frequency range.
When you rectify it, the first phenomenon is that the signal becomes positive. Hence, it is not zero-mean anymore, and thus recovers low-frequencies that were attenuated in the high-pass filtered signal.

A second is the presence of higher frequency terms. Typically, when a sinusoidal component has period $T$, its absolute value doubles its period and generates higher-order harmonics.


Why the carrier frequency does not appear after AM modulation? Okay, this might be a very basic question, but I'm completely new to the field.

I have a $40$ Hz sine wave ($1000$ samples, $fs = 1$ kHz), which I have amplitude modulated to a $4$ Hz sine wave ($1000$ samples, $fs = 1$ kHz). In the time-domain, the result looks like I expected. However, when I calculate the FFT, there are $2$ sidebands but nothing at $40$ Hz anymore. If I increase the FFT resolution with zero-padding I can see something there, but very modest.

Is this expected behaviour? Can someone explain what is happening here and why. Thank you!

EDIT: to clarify, I get two peaks which are at 16 Hz and 24 Hz.

What you did is multiply two sinusoids with each other:

$$z(n)=\sin(n\theta_1)\cdot\sin(n\theta_2)$$

where $\theta_1=2\pi\cdot 40/1000$ and $\theta_2=2\pi\cdot 4/1000$.
The result $z(n)$ can be written as

$$z(n)=\frac{1}{2}\left[\cos(n(\theta_1-\theta_2))-\cos(n(\theta_1+\theta_2))\right]$$

So you get two sinusoids with the sum and the difference of the two original frequencies, i.e. one with $40-4=36Hz$, and the other one with $40+4=44Hz$. This is exactly what you see in your FFT plot. You use a 1000-point FFT, i.e. the frequency at index $i$, ($i=1,2,\ldots,1000)$, is given by

$$f_i=\frac{f_s}{1000}(i-1)=i-1,\quad i=1,\ldots,1000$$

where $f_s=1000Hz$ is the sampling frequency. You get peaks at $i=37$ and $i=45$, which correspond to frequencies $f_{37}=36Hz$ and $f_{45}=44Hz$, exactly as expected.

I think you made a mistake in the mapping of FFT indices to frequencies.

By the way, with standard amplitude modulation you indeed get the carrier in the spectrum because the modulated signal looks like

$$s(n)=[1+mx(n)]\sin(n\theta_0)$$

where $x(n)$ is your message, $m$ is a scaling factor, $\theta_0$ is the carrier frequency, and the term $1$ which is added to the signal results in the carrier wave in addition to the sidebands. You just multiplied two sinusoids without this constant term, that's why you don't see any carrier. This was already correctly pointed out in chirlu's answer.


Do I need an I-Q (in-phase/quadrature) mapper in Simulink with the QAM modulator? Would it be correct if I use the Random Integer Generator block with a certain $M$-ary number (say $M=16$) followed directly by the General QAM Modulator Baseband block?
I saw an example where they use an IQ mapper between the two blocks. Is it necessary?

I assumed the QAM block already does the mapping (inherently) because when I read the document it was specified that the input of the QAM modulator block should be an integer between $0$ and $M-1$, with no mention of the IQ  mapping.

You don't need the mapper. The Simulink QAM modulator lets you specify a vector of constellation points. The input should be an integer (or vector of integers) that select elements from the constellation points.

By default the block specifies a 8-PSK constellation. If you want something different, you need to provide the vector of constellation points.


Is a filter/control transfer function with positive phase "causal"? In control we often use transfer functions with positive phase, i.e., a "lead filter" has transfer function

$$G_c(s) = \frac{\alpha \tau s+1}{\tau s+1}$$

(with $\alpha>1$). Since the zero occurs at a lower frequency than the pole, the phase is positive approximately beteen $\omega = 1/\alpha \tau$ (the zero) and $\omega = 1/\tau$ (the pole), with the maximum phase occuring at the midpoint $\omega = 1/\tau \sqrt{\alpha}$. Typically it's used to increase phase at a particular frequency $\omega_c$ for stability.

But how can a system with positive phase actually exist?
Does the positive phase at certain frequencies mean that if I remove the lead filter from the control loop and input a sinusoid, the output will occur before the input?
How is this related to causality and realizability?

In general, from the positiveness of the phase nothing can be concluded about the causality of the corresponding system. Note that the phase of the given system is

$$\phi(\omega)=\arctan(\alpha\omega\tau)-\arctan(\omega\tau)\tag{1}$$

which is positive for all $\omega>0$ if $\alpha>1$. Yet the system is causal, which can be easily seen by calculating its impulse response:

$$\begin{align}G(s)&=\frac{\alpha s+\frac{1}{\tau}}{s+\frac{1}{\tau}}=\alpha\frac{s+\frac{1}{\tau}-\frac{1}{\tau}}{s+\frac{1}{\tau}}+\frac{1}{\tau}\frac{1}{s+\frac{1}{\tau}}=\alpha+\frac{1-\alpha}{\tau}\frac{1}{s+\frac{1}{\tau}}\\\Longleftrightarrow g(t)&=\alpha\delta(t)+\frac{1-\alpha}{\tau}e^{-t/\tau}u(t)\tag{2}\end{align}$$

where $u(t)$ is the unit step function. Obviously, we have $g(t)=0$ for $t<0$, i.e., the system is causal.

Note that the system's phase is not only positive, but it also has a positive slope for $\omega<1/(\tau\sqrt{\alpha})$. This means that its group delay (the negative derivative of the phase) is negative in that frequency range, as shown in the figure ($\alpha=2$, $\tau=1$):

So even a negative group delay in some frequency range says nothing about the causality of a system.

A rational transfer function with all its poles in the left half-plane (as is the case in your example) and a region of convergence $\text{Re}\{s\}>0$ always corresponds to a causal and stable system, even if its phase is positive or its group delay is negative in a certain frequency range.


How to downscale an image with interpolation by a factor larger than interpolation kernel size? I'm experimenting with image resizing techniques and algorithms. Specifically, I'm significantly downsizing images, e. g. from 2048x1536 to 64x48 - 32 times.
Now, say I'm using a 4x4 kernel. Right now I'm upscaling the image to be by a factor of power of 4 larger than the target (2048x1536 becomes 4096x3072), and then process the image 3 times, each time obtaining a new image 4 times smaller than the source. Seems ineffecient. Is there a trick to it?

It is arguably inefficient. You could achieve the same results with a single processing step using a single kernel. You can do this by convolving all of the kernels together, then applying your resize using this larger kernel.

A similar process to this is used in most digital painting packages. The resizing usually offers some fixed kernel sizes (e.g. bilinear / bicubic - these may be referred to as 'smoother' / 'sharper') as well as some form of intelligent resize (sometimes called 'smart' or 'best quality') where the resampling kernel size is a factor of the image size.

If you need more information you should add some details about what precisely you are trying to achieve.


Texture-like measures for quantifying density of data in binary images Consider the following black-and-white image. It depicts a freehand sketch.

I wish to characterize the "density" of sketch strokes. For e.g. the hair strokes are densely grouped together. So are strokes near the wrists and the necklace stone. Other strokes are somewhat scattered and "far", e.g. the nearly vertical strokes depicting the dress.

Is there a good measure (texture-like ?) which can be used to quantify the above notion ?

One of the most obvious that comes to mind is MPEG-7's Edge Histogram Descriptor.

In this, you divide the image into blocks and consider the edge pixels across different angles forming the complete local histogram. This is organized the same way as a Gabor transform features at angular orientation at different angles. This is fairly a very good indicator to classify and search pictures with similarity in their structures and hence edges.

See this reference 2 which is a great introduction to MPEG 7 color and texture. See Reference 3 for a more detailed coverage for the Local Edge descriptors.

Thomas Sikora, The MPEG-7 Visual Standard for Content
Description—An Overview IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS FOR VIDEO TECHNOLOGY, VOL. 11, NO. 6, JUNE 2001
B. S. Manjunath, Jens-Rainer Ohm, Vinod V. Vasudevan, and
Akio Yamada Color and Texture Descriptors IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS FOR VIDEO TECHNOLOGY, VOL. 11, NO. 6, JUNE 2001
Dong Kwon Park, Yoon Seok Jeon, Chee Sun Won Efficient Use of Local Edge Histogram Descriptor Proceeding MULTIMEDIA '00 Proceedings of the 2000 ACM workshops on Multimedia Pages 51-54


Camera Calibration of a Depth Camera For a regular RGB camera, reading the wikipedia page for camera resectioning/calibration here, they divide the paramters into two categories:
Intrinsic: - focal length, scale factors in x and y, and principal point
Extrinsic: - the rigid body transform of the camera relative to some origin point
My question is how does camera calibration for a depth camera relate to these parameters, for example a Kinect v2 depth camera.
Edit: For example in the SUN RGB-D dataset, there are two files associated with each png image/ png depth image pair. An intrinsic file, and an extrinsic file.
The intrinsic file is a 3 by 3 matrix, example:
529.500000 0.000000 365.000000
0.000000 529.500000 265.000000
0.000000 0.000000 1.000000
From what I understand, using 1-indexing, entry (1,1) and (2,2) according to the camera resectioning article on wikipedia are the focal length multiplied by scale factors relating "pixels to distance" (any elaboration on this appreciated). Whereas entries (1,3), (2,3) represent the coordinates of a principal point. Is this principal point in 2-D or 3D space and what do these numbers mean?
The extrinsics file/s (for some reason there are multiple files sometimes) are
a 4 by 4 matrix, example:
0.979445 0.201365 0.011827 0.000000
-0.201365 0.972640 0.115861 0.000000
0.011827 -0.115861 0.993195 0.000000
Which is a rigid body transform.

Pretty much the same.
The intrinsic calibration tells you the relationship between stuff "out there" relative to the camera's frame of reference and the image the camera makes.  The extrinsic calibration tells you how to get from the camera's frame of reference to whatever your "global" frame of reference may be.


MATLAB Image-Processing Toolbox: asking for a roadmap! I plant trees in the north of BC for a living and have next to
    none sp experience, please be gentle :^]

TLDR:
The problem is essentially this:
 tracking and calculating the "height" of three differently-colored regions whose proportions change relative to each other over a number of images using k-means segmentation and color-histogram vector.

(see fig. below)

It's a little chemistry experiment which kicks off with a transparent beaker full of light-blue-aquamarine foam, then a transparent oil is added and over the next few minutes the light blue foams condenses into a distinctively-darker blue homogenous liquid at the bottom of the beaker and there is no light-blue foam left whatsoever. The boundary between the foam and the emulsifier(oil) is.. well.. foamy at all times as i tried to indicate in the drawing, but i think a blur with an appropriate kernel should solves that pretty well.
 We have the full control over the setup(any lighting, any angle and number of photos, etc.)

They want me to grab the MATLAB image-processing toolbox which I, unfortunately, happened to have and given a set of oh let's say 360 snapshots of the above experiment plot the "height" of the boundary lines between the oil, the light-blue foam and the dark blue liquid at the bottom against time. Without sticking a ruler next to the vessel! Bummer!

How would a signal-processing guru like yourself have gone about this? Terms like "k-means segmentation" and "dominant vectors of a histogram" prevailed in the conversation with the engineer who has imparted this task to me, but knows less MATLAB than I do. He has also made it abundantly clear that he wants it done in MATLAB. Woulda been a simple problem in opencv, but oh well!

Could you possibly throw like a bullet-list of functions/ sequence of steps one would undertake in solving something like this? So i know what to research at least.
 It has to be k-means, histogram and matlab unfortunately.
The difference between going this road and any other for me is like between getting laid and getting shot, folks. I count on your understanding and am indebted in advance!

artie

If the distinct contents/ingredients really have different colors than you just have to take a picture, always in the same location/distance, and then count the amount of pixels for each color.

You can transform your RGB image to the HSV color space to make it easy to identify specific colors.

Some pseudocode:

image_hsv = rgb2hsv(image_rgb)
greenmax = 120+10
greenmin = 120-10
num_green_pixels = sum( and( image_hsv[H]>greenmin, image_hsv[H]<greenmax ) )


How to Extract High Frequency and Low Frequency Component Using Bilateral Filter? What does High frequency and low frequency component signify in an image.
How to separate out the high frequency as well as low frequency component of an image using Bilateral filter.

Similar to one dimensional signals, low frequencies in images mean pixel values that are changing slowly over space, while high frequency content means pixel values that are rapidly changing in space.

For example, the following image has strong low frequency components: You can intuitively see how I simply have a sin-wave propagating at some low frequency.

In contrast, this image below is composed of a sin wave at double the above frequency. Notice how this manifests itself on the image:

Any image can have any number of low and high frequency components together as well. For example, and image like this has both low and high frequency components:

You can see how you have a low-frequency 'trend', but also a lot of high frequency detail across the image. (Roughly speaking, in images, sharp transients like edges correspond to high frequencies, while long unchanging spaces correspond to low frequencies).

So now, the point of an application of a bilateral filter (which is simply a convolution of your image with a gaussian kernel), is to remove high frequency components, and retain your lower frequency components. So in this case, what happens if we convolve the above image, with a gaussian (bilateral) filter that looks like this?

If we use this as our kernel, we will get the following result:

In other words, we have removed the high frequency components of the original image, but kept the low frequency components, since the variance of the gaussian kernel was chosen appropriately.