Can you give example applications using single bin sliding DFT?

I do not want to do FFT since I only need to see a single frequency bin and also FFT is not suited well for real-time applications. I see that there are sliding Goertzel and sliding single bin DFT methods

FFT is just an implementation of the DFT. So, for the same reasons you don't want the FFT you wouldn't want any other DFT.

Goertzel is basically a way of getting the same result any DFT would have given you, but only for a single bin, exactly. Goertzel's latency, when implemented in parallel hardware, is the same as the overall DFT's, just with less ressource usage. In software, you'd use a bunch of operations.

To detect a single tone, my recommendation here would simply be a correlator (which happens to be the same as a matched filter) with the tone of interest. Basically, take the tone you're interested in, and use it as filter taps (apply windowing to taste). You can easily shift in single samples into such a filter. If you don't window (==rectangular window) and convolve a complex sinusoid with a whole number of samples as period, you get what you'd get as a single bin DFT: it's the same formula.

The whole truth, though, is that you don't really get a useful estimate for the presence of a periodic signal with every new sample – no matter how you twist and turn this, it's always a rate-reducing operation ($n$ samples go in, 1 goes out).

So, you might want to consider your definition of "real-time detection". How could one sample change the fact that there's a let's say 200-samples-long periodicity in the signal?

Depending on what you implement this on (an FPGA? ASIC? DSP? MCU? PC-style computer?), you'd want to pick a different implementation, and you'd be getting samples in buffers, anyway – so, I'll go out on a limb here and claim that if you are writing software, you can probably get away with a simple narrowband FIR (and you'd tune that bandwidth and transition width to make your detector robust to slight deviations in frequency), and it's often most efficient to implement that using fast convolution based on the FFT.


How to Calculate Gaussian Kernel for a Small Support Size?

Gaussian Kernel is made by using the Normal Distribution for weighing the surrounding pixel in the process of Convolution.

Since we're dealing with discrete signals and we are limited to finite length of the Gaussian Kernel usually it is created by discretization of the Normal Distribution and truncation.

I created a project in GitHub - Fast Gaussian Blur.
It uses many methods to approximate the Gaussian Blur Filter and evaluate their speed and quality.

The most classic method as I described above is the FIR Truncated Filter:

https://github.com/RoyiAvital/FastGuassianBlur/blob/master/ApplyGaussianBlur.m

Regarding small sizes, well a thumb rule is that the radius of the kernel will be at least 3 times the STD of Kernel.
If you chose $ 3 \times 3 $ kernel it means the radius is $ 1 $ which means it makes sense for STD of $ \frac{1}{3} $ and below.

Then just fill your Kernel By:

gaussianBlurRadius  = ceil(stdToRadiusFactor * gaussianKernelStd);
vGaussianKernel = exp(-([-gaussianBlurRadius:gaussianBlurRadius] .^ 2) / (2 * gaussianKernelStd * gaussianKernelStd));
vGaussianKernel = vGaussianKernel / sum(vGaussianKernel);


Kalman for 3D position and 1D orientation

Reading your link in more detail, it says:

  Acceleration ($\ddot{x}$ & $\ddot{y}$) as well as position ($x$ & $y$) is measured.

which, to me, is not a constant acceleration.  What it is doing is a discretized changing acceleration model, assuming that the acceleration is constant between each sample period.

Now to your question: is inclusion of $\phi, \dot \phi,$ and $\ddot \phi$ in the state vector valid?

The answer is yes, but you need to include the dynamics of how those states interact with the other states.

OK, so now you've added the extra information, I can take a stab at it.

It should work with the following provisos:

The noise distribution on the angle $\phi$ will be uniform not Gaussian.  I can't find a good reference, but this discussion on comp.dsp mentions it.  That means the KF equations for $\phi$ won't be quite right.  I don't think this is a big deal.
There will be some coloration between the noise on $\phi$ and the noises on $x$ and $y$.  I'm not sure what this should be, but it may mean you need to introduce some off-diagonal terms in the $Q$ matrix in the KF formulation.


As long as the two signals are simply delayed versions of each  other the delay is simply given by the point in time where maximum of the cross correlation occurs.

If the signals are uncorrelated, then the question is meaningless.

The tricky part is when the two signals are partially correlated or are filtered versions of the same original signal. In this case you will first need to define what exactly you mean by "delay". One possible way to approach this would be to time shift the signals against each other, subtract them and sum the squares. The shift that results in the minimum energy in the difference would be considered the delay. Something like $d = min[sum((x(n)-y(n-d))^{2})]$


Difference between cross correlation definitions (signal processing / random processes)

If the process is stationary it is ergodic, then time averages are equal to probabilistic averages.


This is actually pretty simple. Take a pure-tone sine wave, say

$$s(t) = \cos\left(\omega t + \phi \right).$$

Now, fix a time $t_0$ and find the phase $\psi(t_0)$ of $s(t_0)$. We get

$$\psi(t_0) = \omega t_0 + \phi$$

Now let's vary the frequency $\omega$ and keep the time fixed. The above expression becomes a linear function in $\omega$! Even though we haven't actually "delayed" the wave.

Imagine it visually. You have a point in time (not 0, that will always stay the same) and a sine wave which is, say, shrinking (frequency increases). As it shrinks, different phases of this sine wave will cross your fixed time point for different frequencies.

This is why phase response of a simple delay system is a sloped line, not a horizontal line. Its derivative, however, will be horizontal and proportional to the amount of delay we introduced in the system, i.e

$$\frac{\partial}{\partial \omega} \left( \omega t_0 + \phi \right) = t_0,$$

Which is precisely the system delay.


I don't think that FM has anything to do with this, unless you are using the term "spectrum inversion" in an unusual context.

The original signal is presumably a real signal centered at 1200 MHz with a bandwidth of 100 MHz.  If you sketch its spectrum conceptually, it will have some content from 1150 to 1250 MHz, and a conjugate-symmetric version of the same thing from -1150 to -1250 MHz.  The conjugate symmetry follows from properties of the Fourier Transform and the fact that the original signal is real.  Call the portion of the spectrum at positive frequencies the positive image, and the portion of the spectrum at negative frequencies the negative image.  If your sampling ends up reversing the positions of the positive and negative images, or if you convert the negative image to complex baseband, then that's spectral inversion.

To find out whether your a sampling scheme causes spectral inversion, recall that the spectrum of the sampled signal repeats at intervals of the sampling rate.  Sketch this operation for the proposed sample rates and look at the position of the positive and negative images in the region from $-f_2/2$ to $+f_s/2$.  Does the spectrum have the original sense, or has it been inverted?


Filter memory refers to the effect of past input on the current input. Infinite impulse response (IIR) filters and finite impulse response (FIR) filters are fundamentally different from the perspective of their memory. The memory of an FIR will essentially be the number of taps, while IIR has theoretically infinite memory (although when they are stable, the memory passes under the numerical accuracy of the system beyond a certain point).

More on this http://gregstanleyandassociates.com/whitepapers/FaultDiagnosis/Filtering/Filter-Memory/filter-memory.htm.


Time-invariance, causality and stability of $h(t)$ of four given systems

You should remember the necessary and sufficient condition for (BIBO)-stability of a system described by an impulse response $h(t)$:

$$\int_{-\infty}^{\infty}|h(t)|dt<\infty\tag{1}$$

If you use condition $(1)$ you will find out that actually only one of the four systems is stable.

The question about time-invariance does not make much sense, because if we assume that the systems are characterized by their impulse response then we already imply that the systems are linear and time-invariant. For linear time-varying systems you would need a two-dimensional impulse response, and non-linear systems cannot be characterized by an impulse response.


A good answer might constitute a few chapters in a good DSP text, making this question overly broad.  So I'll sketch an over-simplified skeleton, where you will need to seek a decent rationale elsewhere.

Pretend your captured signal is just a vector (array of samples).

Move your signal of interest "down" to "baseband" (resulting in either 2 vectors (X,Y) or 1 complex vector (if you know complex arithmetic).
Do this by multiplying by a vector of cosines, oscillating at -8.1 MHz, to get vector X, and by a vector of sines, oscillating at the same (negative!) frequency, to get vector Y.  Low pass filter (to remove audio etc.) these 2 vectors with a cut-off somewhere around 2 or 3 MHz (see specs for your video).

The angle of a 2D vector of X,Y pairs (use atan2() to compute) will rotate with a speed proportional to the modulating signal (e.g. the video signal you want).  The speed of rotation is the first derivative of the angle with respect to time.  Or roughly, subtract successive angles, to get a time vector proportional to dA/dt.  Unwrap/rewrap phase discontinuities as needed.  Low pass filter (maybe somewhere around 2 or 3 MHz) to remove noise.  Scale as needed.  Result should be a vector containing your raw video signal (perhaps at the wrong scale, perhaps noisy, perhaps inverted, perhaps with/without gamma: fix if possible).


How do I create a frequency vs time plot?

The time vs frequency resolution is a well-known problem, and there are indeed approaches to overcome it. For audio signals, some of the commonly used techniques include: parametric methods ; adaptive resolution (analyze with various time/frequency configurations and patch the results together - Wen X. and M. Sandler, "Composite spectrogram using multiple fourier transforms") ; wavelets/decompositions on overcomplete bases ; and use of phase information to extract the precise location of frequency peaks (IFgram).

However, it appears that the graph you have shown does not use some of these techniques ; so I suspect this is not what you might be looking for. There appears to be some "smearing" on the horizontal axis (for example at t=1.2s) and this is a sure sign that the analysis has been done with a high overlap between chunks.

Indeed, the chunk duration and the number of analysis frames per second do not have to be linked to each other if you allow frames to overlap. So if you want to use 40ms long analysis frame, your grid does not have to be:

frame 1: t=0..t=40ms ; frame 2: t=40ms..t=80ms

It could very well be:

frame 1: t=0..t=40ms ; frame 2: t=10ms..t=50ms

This overlapping can give the illusion of a higher temporal resolution without reducing too much the FFT window size. Note that this can only help in accurately locating an event on the time axis - it will not help resolving two events close in time... Just like increasing the FFT size might help with identifying the location of a frequency peak, but not with the resolution of two adjacent frequency peaks.


What is the difference between sampling over $-4\pi$ to $4\pi$ and 0 to $8\pi$

Your spectrum is zero at all frequencies other than at the carrier and the side bands. If the magnitude is zero the phase is undefined. Since you use finite precision math, the magnitude in your calculation will not be exactly zero but just a really really small number. So the phase can be computed but what you are seeing just random numerical noise.
As you change the time interval, the noise will change as well, but the magnitude will still be extremely small.


Coherent Sampling And The Distribution Of Quantization Noise

This is most easily understood by examining the error between the true sinusoid and the quantized sinusoid. Let's start by defining a simple sinusoid. For simplicity, I will assume a cosine, but the same arguments apply for a sinusoid with an arbitrary phase offset. The true sinusoid is defined in continuous-time as

$$ x(t) = \cos\left(2 \pi f_{in} t\right). $$

Now, we need to sample this signal with a period of $T_s = \frac{1}{f_s}$ to obtain a discrete-valued (but unquantized) signal

$$ x_D(n T_s) = \cos\left(\frac{2 \pi f_{in} n}{f_s}\right).$$

Using the relationship that you provided in your question, we can rewrite this as

$$ x_D(n T_s) = \cos\left(\frac{2 \pi N_c n}{N}\right). $$

From this it should be clear that this signal lies in two bins, specifically the bins at $N_c$ and $N - N_c$ of an $N$ point FFT. However, this signal is not realizable with finite precision numbers. So, we need to quantize it. Let us assume that we are quantizing to $q+1$ bits using two's complement. A reasonable approximation to the quantized signal is given by

$$ x_Q(n T_s) = 2^{-q} \left\lfloor 2^{q} \cos\left(\frac{2 \pi N_c n}{N}\right) \right\rfloor.$$

Another way to represent this would be to represent the sinusoid added with an error term (i.e., the quantization noise $\epsilon_q(n T_s)$). This is captured in the following:

$$ x_{D}(n T_s) = x_Q(n T_s) + \epsilon_q(n T_s). $$

Now, let's consider what $\epsilon_q(n T_s)$ looks like for the two cases you mentioned. That is, the coprime case and the common factor case. Let's rearrange this slightly so that we can see things a little better,

$$ \epsilon_q(n T_s) = \cos\left(\frac{2 \pi N_c n}{N}\right) - 2^{-q} \left\lfloor 2^{q} \cos\left(\frac{2 \pi N_c n}{N}\right) \right\rfloor. $$

Now, consider the case where $N_c$ and $N$ share a common factor, $M$. Let's say $M$ is the greatest common factor. In this case, $N_c = M K_c$ and $N = M K$. Using this and the expression derived for $\epsilon_q(n T_s)$, we can show that

$$ \epsilon_q(n T_s) = \cos\left(\frac{2 \pi K_c n}{K}\right) - 2^{-q} \left\lfloor 2^{q} \cos\left(\frac{2 \pi K_c n}{K}\right) \right\rfloor, $$

and from this, we can see that

$$ \epsilon_q(n T_s) = \epsilon_q(n T_s + K T_s) $$

since

$$ \cos\left(\frac{2 \pi K_c n}{K}\right) = \cos\left(\frac{2 \pi K_c (n + K)}{K}\right).$$

In other words, the sequence is periodic with period $K < N$ since $K$ and $K_c$ are coprime. This means that all of the information in this signal can be captured with a $K$ point FFT. This, in turn, means that all of the energy from this signal is contained at the equivalent bins in the $N$ point FFT. These are the bins at $0, M, 2M, ... (K-1)M$. All of the other bins are zero-valued. This is why you see a discrete spectrum as shown in the slides you mentioned.

Now, if $N$ and $N_c$ are coprime, then the sequence is non-sinusoidal and aperiodic (or viewed another way, it is periodic with period $N$). There is no shorter period embedded in the sequence as there was with the common factor case. Therefore, all bins are required to represent the signal.


The plot is of $$\mid X\left(i\omega\right) \mid = \sqrt{\left(\frac{1}{a+j\omega}\right)\left(\frac{1}{a-j\omega}\right)} = \frac{1}{\sqrt{a^2 + \omega^2}}$$

against $\omega$

In particular $\omega$ can be equal to $-a$.

This checks out with Wolfram alpha


Is the value of a frequency bin of a DFT-output the average of the 'real' frequency values within that bin's range?

I'll start with a counter-example:

Assume your continuous-time signal is a sinusoid at freqeuncy $f_0$. You then sample that signal, and use a rectangular window to select $N$ samples. Feed those $N$ samples to your DFT, and you will see non-zero values in the whole spectrum (unless in the special cases where $f_0$ is a multiple of $f_s/N$).

So you start with a signal whose energy is concentrated in a single frequency, and you end up with non-zero values in all DFT bins.

To understand what is happening:

Start with a continuous-time signal $x_c(t)$.
(Optionally) filter the signal to make it band-limited to $f_s/2$.
Sample the signal to get $x[n] = x_c(n/f_s)$.
At this point, the spectrum is given by the sampling theorem. You will have spectrum aliasing if the sampled signal was not band-limited in step 2.
Now apply a window of size M. The window type depends on your application (Rectangular (ugly), Hann, Hamming, Blackman, Flat top, etc.).
At this point, your spectrum is that of the sampled signal, but convolved with the window (depending on window size and type, you will get some energy spread to other frequencies). This spectrum is $X_w(e^{j\theta})$.
The DFT of the windowed signal will return samples of its spectrum, at frequencies $Y[k] = X_w( e^{j 2\pi k/N})$.


How to compute generator matrix from a parity check matrix?

With forward-error-correcting coding, one is working in a finite field, typically the field of two elements denoted by GF$(2)$ or $\mathbb F_2$. So, there are no fractional numbers
and no fancy methods such as singular value decomposition: you use bit-by-bit XOR additions of the rows of $H$ and Gauss-Jordan elimination to reduce $H$ to row-echelon form $[P_{(n-k)\times k} \mid I_{(n-k)\times(n-k)}]$.

Then, set $G = [I_{k\times k} \mid (P^T)_{k\times(n-k)}]$ and you are done. (For
nonbinary fields, use $[I \mid -P^T]$).  Note
that all arithmetic in the verification $HG^T = 0$ is also finite field arithmetic
with $1\cdot 1 =1$ and $1+1 = 0 = 1\oplus 1$ for the case of GF$(2)$.


FYI... STASM uses OpenCV's face detector (albeit very slightly modified for better image-boundary detections). If you really want a better face detector then I recommend hunting down detectors from the FDDB dataset challenge (most probably proprietary) or create your own custom detector to fit your needs.

If you still plan on going the route of verification via eyes, then OpenCV has a cascade for that, see http://docs.opencv.org/doc/tutorials/objdetect/cascade_classifier/cascade_classifier.html


time scaling and shifting of cosine in Fourier transform

Using the shifting and scaling properties of the Fourier transform is a rather complicated way of computing the Fourier transform of $\cos(at+b)$. A more straightforward way is to realize that

$$\cos(at+b)=\frac12\left[e^{jat}e^{jb}+e^{-jat}e^{-jb}\right]\tag{1}$$

From the Fourier transform pair

$$e^{jat}\Longleftrightarrow 2\pi\delta(\omega-a)\tag{2}$$

we immediately obtain

$$\cos(at+b)\Longleftrightarrow \pi\left[e^{jb}\delta(\omega-a)+e^{-jb}\delta(\omega+a)\right]\tag{3}$$

But having said that, your result is correct. You just need to simplify it appropriately. In order to do so, you need to know two things. First,

$$\delta\left(\frac{\omega}{a}\right)=|a|\delta(\omega)\tag{4}$$

and, second,

$$f(\omega)\delta(\omega-\omega_0)=f(\omega_0)\delta(\omega-\omega_0)\tag{5}$$

for any $f(\omega)$ that is continuous at $\omega_0$.

From $(4)$ you get

$$\delta(\omega/a-1)=|a|\delta(\omega-a)\tag{6}$$

And from $(5)$ you get

$$e^{j\omega b/a}\delta(\omega-a)=e^{jb}\delta(\omega-a)\tag{7}$$

Applying similar transformations to the second term in your result leads to the form given in $(3)$.


Segmentation of small artifacts - Image Processing

If this is really what your input looks like then you should be able to handle this with the standard morphological image processing operations.

Assuming foreground as white and background as black, it basically comes down to a morphological open with a big enough structure element to elide the tails. The only trouble is that doing so will not work given that the input has holes so large the lines in the body of the shape itself would also get nuked by the open operation. One way to handle this is to perform the open on the figure without holes and then put the holes back in.

For example, in C++ using OpenCV:

void Test()
{

    Mat src_bw;
    cv::threshold(
        imread("some_file.png", CV_LOAD_IMAGE_GRAYSCALE),
        src_bw, 0, 255,
        CV_THRESH_BINARY | CV_THRESH_OTSU
    );

    imshow("src_bw", src_bw);

    // perform contour extraction using RETR_CCOMP i.e. we are assuming input doesn't have nested islands/holes just simple holes.
    std::vector<std::vector<Point> > contours;
    std::vector<Vec4i> hierarchy;
    findContours(src_bw.clone(), contours, hierarchy, RETR_CCOMP, CV_CHAIN_APPROX_SIMPLE, Point(0, 0));

    // iterate over the contours, draw the top-level contours into "blobs"
    Mat blobs = Mat::zeros(src_bw.size(), CV_8U);
    for (int i = 0; i < contours.size(); i++)
        if (hierarchy[i][3] == -1)
            drawContours(blobs, contours, i, 255, CV_FILLED, 8, hierarchy, 0, Point());

    imshow("blobs", blobs);

    // do an 'open' on blobs which will remove the artifacts.
    auto shape = MORPH_ELLIPSE;
    Mat strel = getStructuringElement(MORPH_ELLIPSE, Size(5,5));
    Mat blobs_opened;
    morphologyEx(blobs, blobs_opened, MORPH_OPEN, strel);
    imshow("blobs_opened", blobs_opened);

    // draw the holes back into the output
    Mat output = blobs_opened.clone();
    for (int i = 0; i < contours.size(); i++)
        if (hierarchy[i][3] != -1) {
            drawContours(output, contours, i, 0, CV_FILLED, 8, hierarchy, 0, Point());
        }

    imshow("output", output);
}

yielding the following:


How to determine the window size and weights in Weighted Moving Average (WMA), given desired cut-off frequency?

Your normalized window is given by

$$w(n)=\frac{2}{N+1}\frac{N-n}{N},\quad n=0,1,\ldots,N-1$$

The window satisfies

$$\sum_{n=0}^{N-1}w(n)=1$$

which means that the gain of the corresponding moving average filter is 1 at DC.

For determining the cut-off frequency, we need to compute the frequency response of the window:

$$W(e^{j\theta})=\sum_{n=0}^{N-1}w(n)e^{-jn\theta}$$

After some algebra you get

$$W(e^{j\theta})=\frac{2}{N+1}\frac{1-\frac{N+1}{N}e^{-j\theta}+\frac{1}{N}e^{-j(N+1)\theta}}{(1-e^{-j\theta})^2}\tag{1}$$

Now you need to find the value of $N$ for which the magnitude of (1) at the cut-off frequency $\theta_c=2\pi\frac{3Hz}{50Hz}$ becomes $1/\sqrt{2}$ (-3dB). Since $N$ must be integer you cannot get any desired cut-off frequency, but the given cut-off is approximately achieved by $N=9$, for which $|W(e^{j\theta_c})|=0.698$ (-3.13dB).


But I really wanted to write coefficients and formulas!  Oh well, ok...in simplest terms a low pass filter is a weighted moving average. Consider charts of stock market prices and look at the minute by minute pricing compared to a "30 day moving average" and ask yourself which plot contains higher frequencies?

By adjusting the weights of each delayed copy of our signal when we combine them in the process of determining the average, we can change the characteristics of the filter in how it attenuates versus frequency.  So consider the filter as a moving window that we can average over (the memory of the filter) and higher frequency components that have many cycles within the window duration will tend to average out over this window length, while lower frequency components will not go fully plus and minus over the duration of our window so will not average out and will pass through.

For all filters analog and digital, to attenuate faster, more delay is necessary. It is amazing how much we can do from just delays and adds with weighting of each of those delayed outputs!

Update: It appears the OP is asking about a high pass filter, while the above description pertains to a low pass filter.  Note that we can take a linear phase (constant delay) low pass filter and subtract it from a delayed and gain (DC gain) matched version of the input (matching the delay of the filter but passing all frequencies equally), we would achieve the desired high pass response.  That may help to see how the same factors related to delay and filter performance apply. (This is not how I would typically implement a high pass filter as we can modify the coefficients directly, but gives insight into the process).


Given a -3dB/ octave filter that makes Pink Noise, how can I make a +3dB/ octave filter that makes Blue/ Azure Noise?

You need to invert the filter, i.e. flip the poles and zeros.
The implementation that you reference is fairly awkward and will require a decent amount of math work to invert: you need to write the Z-transform for each first order section, add all the fractions into a single fraction and calculate the zeros of numerator polynomial.
An easier way would be to use this implementation https://ccrma.stanford.edu/~jos/sasp/Example_Synthesis_1_F_Noise.html and simply swap $A$ and $B$


Calculating 1/3 Octave Spectrum from FFT / DFT

Dr. Manuel Kuehner,
  You are close. You need to take the square root of the linear values squared.

$$P_{\mbox{total_linear}}=\sqrt{p_1^2+p_2^2+...}$$

$$P_{\mbox{total_dB}}=20 log_{10}\left( P_{\mbox{total_linear/20E-6}} \right)$$

FYI: I wrote a MATLAB function to do exactly as you request. It is here

Looking to read? See page 16 of this book:
https://www.amazon.com/Noise-Control-Engineers-Harold-Lord/dp/0070387389


how can I plot a power spectrum for 150 range bins with 512 fft points

By default, matlab's plot command clears the plot area before drawing the new curve.

To prevent this, and allow multiple curves to appear on the same figure, you can activate the hold mode:

hold off; % 1st plot will clear figure
for i = 1:length(data header.nrgb)
   plot(freq,arraydata(i,:)+height(I));
   hold on; % set hold mode for subsequent figures
end

Alternatively you can stack multiple copies of freq and the plot command will draw one curve per column. Since your data is organized by rows, you'll need to transpose your data matrix:

% just make sure freq & height are column vectors. This may not be necessary if they already are.
if size(freq,1) < size(freq,2)
   freq = transpose(freq);
end
if size(height,1) < size(height,2)
   height = transpose (height);
end

% plot one curve per column on same figure
N = length(data header.nrgb);
plot(repmat(freq,1,N), transpose (arraydata)+repmat(height,1,N))


Why does the flat part of the BER curve appear after the Viterbi decoder (MatLab)?

while numErrsSoft < 100 && numBits < 1e7

Nope, you can't stop after 10⁷ bits if you're after error rates < 10⁻⁶. That means you expect to see less than 10 bit errors. That's not a sufficient statistic!

So, this is where I know your curve becomes unreliable, anyway.

You might still be running into some error flow due to the implementation, but I can't look inside Matlab's implementation of neither the QAM decision nor its Viterbi, so these are "meh" at best for scientific considerations about the nature of Viterbi.


Signal matching: How to quantify the difference between two signals of varying frequency

For a relatively easy answer you can use a short-time Fourier transform and use the largest value as the estimate of your sinusoid frequency.  The basic idea is that you do short Fourier transforms at different points in time to see how it changes.  The problem with that is that you have a tradeoff between time resolution and frequency resolution.  Shortening the Fourier transform improves the time resolution at the expense of the frequency resolution.  Conversely, making the Fourier transform longer improves the frequency resolution at the expense of the time resolution.

The frequency resolution can be improved by using quadratic interpolation to get a finer estimate of the true sinusoid peak.

Another way to go that can get you very good frequency estimation at every point in time is with a phase-locked loop.  Phase-locked loops are a little trickier to set up though, because you have to get the feedback gains right.  There is also a tradeoff in regards to responsiveness (akin to time resolution) and resistance to noise.


How to interpret these different Fourier analysis of this audio signal?

If you look at the Rectangular window the best its rejection gets is about 40 dB. So that behavior, especially obvious in the bottom plot, for the rectangular window is to be expected.

I don't know for sure if this explains everything but look at the level between your peak signal and the high-frequency components. There is almost 60 dB of rejection there. I've always heard that a good rule of thumb is to get 60 dB of rejection from your filters. I know that's not the full 96 dB offered by 16-bit music, but I bet you'd be hard pressed to actually hear that. Of course, there's the silliness of having music at that sample rate. Humans just can't hear anything above around 20 kHz, give or take. This article gives a good summary of the issues. It also brings up a good point, that that sample rate has the potential to pick up harmonics and other high-frequency distortion caused by equipment and electronics, despite the fact that we can't hear it. Perhaps something like that is going on?


How to model a generic low frequency signal?

Note that you can also use Parseval's theorem. Assuming that $f(t)$ and $g(t)$ are both real-valued we have:

$$\int_{-\infty}^{\infty}f(t)g(t)dt=\frac{1}{2\pi}\int_{-\infty}^{\infty}F(\omega)G^*(\omega)d\omega\tag{1}$$

If $g(t)$ is a narrow-band low-pass signal, i.e., if we can model $G(\omega)$ as $G(\omega)\approx G(0)\text{rect}(\omega/\Delta\omega)$, and if we can assume that $F(\omega)$ is approximately constant (equal to $F(0)$) in the narrow interval $\Delta\omega$, then we can write

$$\frac{1}{2\pi}\int_{-\infty}^{\infty}F(\omega)G^*(\omega)d\omega\approx\frac{1}{2\pi}F(0)G^*(0)\Delta\omega\tag{2}$$

where $\Delta\omega$ is twice the bandwidth of $g(t)$ (because the bandwidth is defined as the support at positive frequencies). With

$$F(0)=\int_{-\infty}^{\infty}f(t)dt$$

and

$$G^*(0)=G(0)=\int_{-\infty}^{\infty}g(t)dt$$

we obtain from $(1)$ and $(2)$

$$\int_{-\infty}^{\infty}f(t)g(t)dt\approx\frac{\Delta\omega}{2\pi}\int_{-\infty}^{\infty}f(t)dt\int_{-\infty}^{\infty}g(t)dt\tag{3}$$


In the first channel, you can find $\alpha>0$ by first transmitting a value $x>0$ known by the receiver; then, the receiver can find $\alpha = r/x$. For example, if $x=0.5$, and $r$ is received, then $\alpha = r/0.5$. Once you know $\alpha$, the second symbol allows you to transmit at an infinite rate. Say you want to transmit the bit sequence $b_1\,b_2\,b_3\ldots$. Let $x=0.b_1\,b_2\,b_3\ldots$ in decimal, and you're done. Note that the rate is infinite because you're transmitting an infinite number of bits in two milliseconds: $$\text{Rate} = \frac{\text{information}}{\text{time}} = \infty.$$

For the second channel, note that you can't use that capacity formula, because it assumes that $n_i$ is Gaussian, which is not your case. Rather, you can transmit with zero error probability by defining a set of amplitudes for $x$ such that the noise cannot introduce an error. As a starting point, consider this scenario. You let $x \in \lbrace 0.25, 0.75 \rbrace$, and define a mapping such that a bit zero is transmitted as $x = 0.25$ and a bit one is transmitted as $x = 0.75$. This system has zero errors, since the receiver can apply a decoding rule with threshold $0.5$. It is not possible for the noise to confuse the receiver in this case.

So, what you need to do is figure out how many different values can $x$ take such that the noise cannot introduce errors. That will lead you to the capacity in bits per symbol. Since you know the symbol rate, you can also find the capacity in bits per second.


digital IIR filter delay and oversampling

Note that a given digital filter impulse response $h[n]$ will have a corresponding equivalent analog frequency response $H(e^{j\omega})$ through the sampling relations; hence with the given sampling rate Fs.

This means that if you change the sampling rate of the input of this filter, then the effective analog filter response will also be changed, probably that's not what you want to have. So you have to redesign the digital filter to recreate the desired analog filter response.

In particular, when you increase the sampling rate by two, then the digital filter $h[n]$ should be redesigned such that its digital frequency response cutoff frequency is halved, so that the associated analog cutoff frequency would remain the same.

That means that, since the digital cutoff frequency is halved, the filter impulse response is doubled in length, the delay is doubled too. But don't worry, as the sampling period is also halved, the associated analog delay remains the same.

Of course, delay here refers to that associated with the decay time of the digital IIR impulse response $h[n]$.


Processing and feature extraction from very short signals

Though problem! However, signal processing might have a tool at hand: It's called compressive sensing, and reduces the number of samples you need to sub-Nyquist-rate levels.

It's a bit nonsensical to derive the math behind compressive sensing for signals that are a sum of sinusoids (that being the only thing that I could do without literature) if that doesn't describe your signals overly well. However, the idea is:

If you want to know about everything happening with less than a given frequency over a specific time, you'll need all of the samples at a constant rate (typically, twice the highest frequency for real-valued samples). Then you can, for example, do a base transform to get the samples from time domain to frequency domain using the DFT.

For the following, assume that you could, if you had the "full, long, evenly-sampled" signal $\mathbf S \in \mathbb R^N$, define a base transform $\mathbf T$ to transform the time samples vector to a result vector giving you the information you want $\mathbf Y = \mathbb R^N$. That base transform would then simply be a matrix:

$$ \mathbf Y = \mathbf S^T \mathbf T$$

If what you want is sparse, ie. the dimensionality of the result vector would be much smaller than the observed samples vector, the transform matrix wouldn't be a square full-rank base transform, but a matrix that would have $N$ rows, but fewer columns.

The point here is that most of the values in $\mathbf T$ would be pretty small -- setting those to zero will incur an error, that "feels" like noise, but given a nonzero upper limit for allowable distortion, you can find a sparse matrix that still can allow for reliable conclusions about $\mathbf Y$.


Can you decimate / downsample a signal in frequency domain just like you can interpolate / upsample it?

Your interpolation operation is an upsampling operation with subsequent low-pass filtering. This translates to zero-padding in the frequency domain. Given your original signal $x[n]$, you upsample it to become

$$y[n] = \begin{cases}x[n/M]& n/M \in \mathcal{Z} \\ 0 & \text{else}\end{cases}$$

to become an M-times upsampled signal. Then, your interpolated signal is

$i[n] = h[n] * y[n]$ where $h[n]$ is a Dirichlet kernel/periodic sinc function with zeros of distance $M$. In the spectral domain, this corresponds to

$$Y(f) = \sum_{m=0}^{M-1}X(M(f-m))=X(Mf)$$
I.e. a M-times periodic repetition of its spectrum (Note that a periodic spectrum creates the discrte time-domain signal).  Here, I assume $f=0...1$ for the discrete case and the last equality follows from periodicity of $X(f)$ here. Also note that your original spectrum is narrowed, since your sampling frequency increased.

Now, you filter out the spectrum repetitions by multiplying with a rect-function in the frequency domain. This corresponds to convolution in time with the sinc function.

Now, to do the reverse (i.e. decimation), this is what happens: Let $x[n]$ be your original signal you want to decimate. You calculate

$y[n] = x[Mn]$ and this corresponds to the following in the spectrum:

$$Y(f) = \sum_{m=0}^{M-1} X((f-m)/M)$$

I.e. your original spectrum is repeated M times overlapping each other and summed together. This in general creates aliasing. Hence, a simple removal of higher frequencies in the DFT domain is not the same as taking every Mth sample from the time domain signal. Instead, you need to take the spectrum, sum it up repeatedly (such that it becomes M periods within your full frequency band) and finally perform the ifft of one period of the spectrum.

Your approach works, if your signal is bandlimited enough, if it does not create aliasing, i.e. the bandwidth of the signal is smaller than 1/M. Then, the shifted spectra dont overlap and you can just remove do ifft of the smaller portion. Alternatively, you can see your approach as filtering with an aliasing filter (i.e. lowpass with cutoff 1/M) before downsampling.


Your question is overly broad and fundamentally impossible to answer because one needs preparation and no one knows how prepared you are.

Matlab and signal processing are 2 different things. It is used a lot by people who do signal processing but they are not the same.  Learning Matlab if you have some programming background and some linear algebra is not hard.  I’ve seen 60 year old die hard FORTRAN programming engineers pick it up in 2 weeks.  They understood the numerics and it was just a matter of reexpressing their knowledge.

Using a tool you don’t know how to use on a problem you don’t understand is very difficult. Did I do it right? Is it the code? Is it the algorithm?  Often it will be both.

FDTOOL is considered easy to use.

The Understanding Digital Signal Processing by Lyons is recommended.

Look through the old posts on this sight, particularly those that are “frequent”.

Some people learn in a line. Other people learn in a circle. They revisit concepts after learning new ones. Circles have no start and end. I think everyone is a bit of both.


Watermarking and steganography are very related to each other. While steganography is interested in capacity and invisibility, watermarking is to provide robustness to protect hidden data against any attacking methods like compressing, cropping, noise addition etc. LSB Coding, Phase Coding, Spread Spectrum, Echo Hiding are basic methods you should search about.

I suggest you to read these articles for the basics:

Audio Watermarking Techniques, Hyoung Joong Kim
Robust and Reliable Audio Watermarking Based on Dynamic Phase Coding and Error Control Coding

First article is very helpful to start. I advice you to try Spread Spectrum method, it was very robust against mp3 compression as I have coded myself.  I had answered my own question here trying to explain steps of Direct Sequence Spread Spectrum Method to embed data inside audio files. It may also help you.


Maintain constant total amplitude of sum of sines

Hopefully this will be some help.

As you have stated it, this is a difficult problem to formulate and thus get an analytical solution.  First, are you talking continuous or discrete?  Makes a big difference.

Anyway, you want to formulate it as an optimization problem.  Something like "minimizing the RMS" can be expressed as an evaluation function.  "All the peaks having the same amplitude" means have to measure the value of the function with the condition applied to the derivative (=0).  So, your evaluation function (a function of your free variables) could be something like:

$$ F(\phi_0,\phi_1,\phi_2,...,\phi_{N-1}) = \sum_{p=1}^{P} \left( f(t_p) - \bar f \right)^2 $$

Where $f(t)$ is your function, evaluated when $f'(t)=0$, and $\bar f$ is the average height of all those peaks.

This evaluation function is going to be periodic in all dimensions.  There are going to be lots of local minima.  Your only recourse (as I can see it) is to do a coarse sweep of your whole domain, then a local optimization search on the best candidates.


First of all, there is something wrong in your computation of the sine wave. It should be something like:

volume * sin(2 * pi * frequency * n / SAMPLE_RATE);

At the moment, your code is generating very high frequencies and what you are hearing are their mirror images.

This is the very first thing to fix. Other problems in your code include:

Resetting the variable n every time the frequency is changed - this creates discontinuities ("clicks") in the waveform.
Resetting the variable n every time it reaches the sample rate. You are not sure the waveform has completed a series of full cycles because the frequency might not be an integer.
Directly computing the phase by multiplication rather than by integration of the frequency. This causes problems when frequency is modulated.
Computing the sin function directly, which is an expensive operation on embedded hardware, that will have to be replaced by a table lookup + interpolation.

Also, make sure that your samples are scaled appropriately with respect to the scale of your DAC (for example, if this is sent to a 12-bit DAC with unsigned codes, volume must be lesser than 2047, and you should add a 2048 offset).

If you get it working, you will have a pure tone - I wouldn't call this tone "unpleasant", but it is static and not reminiscent of any music instrument.

There are many sound synthesis techniques to try, an overview of which wouldn't fit in this post. In short:

Subtractive synthesis works by generating harmonically rich tones with simple waveforms (sawtooth, square), and then shaping them through filters. It is not really good at producing "realistic" sounds, but it is relatively easy to implement (the trickiest step is bandlimited waveform generation, read about minBLEP), the sounds it produce are expressive, and it is easy to understand how each synthesis parameter relates to each aspect of sound. Used by synthesizers in the 65-85, with a revival since 95 due to the popularity of the sounds they produce in electronic/dance music. A way of "cheating" to get subtractive synthesis sounds more realistic is to use single cycle waveforms of real instruments in place of oscillators, a technique known as wavetable synthesis.
Additive synthesis directly builds up a complex tone by summing sinusoidal partials at various frequencies/amplitudes. Given the number of partials involved and the complexity of their variations over time, it is a rather unpractical method.
FM synthesis works by interconnecting simple sine wave oscillators which modulate each other's frequency. It is surprisingly good at imitating natural sounds (esp. brass, metallophones), but it is very difficult to understand how its parameters (frequency ratios and modulation indices) relate to meaningful musical parameters (sound "brightness" for example). Had a golden era in the eighties and faded away.
Sampling uses recordings of real instruments, and adjusts the playback speed of the recording (aka sample) to emulate the different notes. With this method it is trivial to get something that is immediately recognizable as a piano or guitar or coughing Ferris Bueller - though if you wish to capture all the playing nuances of an instrument you will need an increasingly large number of recordings. Currently the dominant form of sound synthesis in many applications, driven by the lowest cost of data storage.
Physical modeling attempts to simulate the physical processes happening in musical instruments - using techniques such as finite element models or digital waveguides. This process can be computationally expensive and/or require a complex modeling and data analysis effort. The results can be very expressive in reproducing all playing nuances of an instrument. Currently quite rare, except for some very narrow application (Piano synthesis).

It is really a case of "realism/expressivity/simplicity, pick any two".

I am not sure about your goals - how many instrument sounds should your system provide? do you also want to modulate other parameters of sound besides pitch and volume with the accelerometer data? I am not sure either about your target platform (the ST MCUs go from 24 MHz/32kb flash/no FPU to 168 MHz with FPU and DSP instructions and 1M of onboard flash). Another factor to take into account is polyphony - do you want your system to be able to play several notes at the same time?

I would recommend:

Sampling if you have at least 128kb of flash and need only one sound/instrument, or if you have very limited CPU.
Subtractive synthesis if you have less flash and/or need more sounds/instruments ; or if you need to use external sensor data to easily modulate sound parameters besides pitch/volume.


The flaw in the problem posed starts right from the block diagram  that you have posted: it is not the case in DPSK receiver that the RF signal is multiplied with a delayed version of itself and then passed through a low-pass filter followed by slicing etc.

Instead, the RF signal is "demodulated" through a variant of a standard coherent QPSK receiver whose local oscillator is assumed to be synchronized in frequency to the incoming carrier signal but not necessarily synchronized in phase. Thus, the incoming $$r(t) = g(t)\cos(2\pi f_c t+\theta) + n(t), 0 \leq t < T,$$ where $g(t)$ is the baseband pulse, $T$ the bit duration, $\theta$ the unknown RF phase, and $n(t)$ the AWGN (or band limited AWGN if you like to think that the incoming signal has been bandpass-filtered, e.g. as in a tuned amplifier in an IF stage of a superheterodyne system), is multiplied by $2\cos(2\pi f_c t)$ in the I branch and separately by $-2\sin(2\pi f_c t)$ in the Q branch. Since
\begin{align}
2\cos A \cos B &= \cos(A+B) + \cos(A-B)\\
2\sin A \cos B &= \sin(A+B) + \sin(A-B)
\end{align}
we get that the signal parts of the two mixer outputs are
\begin{align}
2\cos(2\pi f_c t) g(t)\cos(2\pi f_c t+\theta) &= g(t)\cos(4\pi f_c t+\theta) + g(t)\cos(-\theta)\\
-2\sin(2\pi f_c t) g(t) \cos(2\pi f_c t+\theta)&= -g(t)\sin(4\pi f_c t+\theta) - g(t)\sin(-\theta).
\end{align}
The double-frequency terms in the mixer outputs can be eliminated explicitly by low-pass filtering the results before doing anything else, or they get eliminated implicitly when we do matched filtering on the two mixer outputs. Now, matched filtering is typically done by via correlation
(multiply by $g(t)$ and integrate from $0$ to $T$) and thus we effectively get a complex-valued decision variable $Xe^{j\theta} + \mathcal CN(0,\sigma^2)$ at the end of the integration period. We can't make a decision just on this decision variable because we don't know the value of $\theta$; depending on what $\theta$ is, the decision might be wholly bass ackwards from what it should be! But what we can do is save the decision variable and use it $T$ seconds later to answer the question

  Did the incoming signal change phase from the unknown $\theta$ during $[0,T)$ to (also unknown) $\theta+\pi$ during $[T,2T)$, or did the phase stay the same unknown $\theta$ during $[T,2T)$?

The final expression for the error probability of DBPSK is quite simple:
$$P_{\text{e, DBPSK}} = \frac 12 e^{-E_b/N_0}$$
but, as you say, is harder to derive than the $Q
\left(\sqrt{\frac{2E_b}{N_0}}\right)$ error probability for coherent BPSK.

Turning to the question of why the trick of using the coherent BPSK error formula as a guide and simply replacing the SNR for BPSK with the SNR for DBPSK in the error probability formula doesn't give the right answer for the DBPSK error probability, the issue is that the two systems make decisions differently and there is no obvious reason why the error probabilities of the two systems should be given by the same formula in terms of SNR. Note, for example, that $P_e = Q(\sqrt{2E_b/N_0})$ for coherent BPSK but only $Q(\sqrt{E_b/N_0})$ for coherent FSK.

Finally, note that for $x > 0$,
$Q(x)  < \frac 12 \exp\left(-\frac{x^2}{2}\right)$
and so, setting $x=\sqrt{2E_b/N_0}$, we get that
$$P_{\text{e, BPSK}} = Q(\sqrt{2E_b/N_0}) < P_{\text{e, DBPSK}} = \frac 12 e^{-E_b/N_0}$$


Algorithms for detection of fingertips or nails

I would consider using neural network or SVM to fit the model. The difficulty with this approach is that you must collect a lot of data - both positive and negative examples. But you can generate a lot of artificial data(by scaling or rotating images you have already collected). After you collect data you can use "moving window" of few sizes to detect nails on images not present in training set.
I don't know how much time you want to spend doing the project. But implementing efficient learning algorithm yourself is difficult. However there is already library for SVM which I would use for this purpose.


Second order sections and normalized factor

If I understand your question correctly, you want to split the total gain between two or more stages. In this case you just have to satisfy $G=g_1g_2g_3$, where $G$ is the total gain returned by tf2sos() and the $g_i$ are the individual gains per stage.


I think you're just misjudging what the [CHANNEL] block does: this is just saying "all the bits that fall out of the different quantizers/encoders get transmitted".
It's especially not an adder! Most books actually show this just as a "frame packer" or "interleaver".


The resolution of the Blackman-Tukey periodogram estimator is $\mathcal{O}(1/M)$, whereas its variance is on the $\mathcal{O}(M/N)$. This compromise between resolution and variance is to be considered when choosing the window’s length.

I have paraphrased Peter Stoica's beautiful explanation from here (check Section 2.5)


Distribution of energy in frequency bands

Frequency bins corresponding to the frequency bands could be found using the fact that each bin will represent band of $(fs/\mathbf{fft}_{pt})$.
$\mathbf{fft}_{pt}$ is the number of fft points you used for stft.

Hence you can find bins corresponding to your band of requirement.


pole/zero locations for real and imaginary signal

I think there is a confusion between signals and systems here. A signal doesn't have any poles and zeros. Only systems have poles and zeros.

A system is a thing that turns one signal into another one. For example a loudspeaker turns a voltage signal into a sound signal. The voltage going in is a signal, the sound coming out is a signal, but the loudspeaker itself is a system.

Linear Time Invariant Systems can be characterized by their transfer function or their impulse response. The transfer function in the Z domain can be written as a rationale function, which in turn can be characterized through poles and zeros.

If the system has only a real valued input and only a real valued output, the impulse response must be real as well and that in turn will require the poles and zeros to be either real or be conjugate complex pairs.


Example of an LTI system with complex impulse response

Complex band pass filters are used to get a band and its quadrature in one efficiently computed step. The most simple such design is a single-pole resonant bandpass with the discrete transfer function

$$H(z)=\frac{(1-r)\exp(i \phi)}{z-r \exp(i \phi)}$$

where $r<1$ determines the bandwidth and $\phi \in [0,2\pi[$ the band center frequency. If bandwidth and center frequency are chosen to avoid significant negative frequency amplitudes, the resulting signal is practically analytic.

Strongly bandlimited and analytic signals have a complex magnitude envelope that is well behaved and can be interpreted as the instantaneous amplitude of the band. They also have a well defined phase derivative that can be interpreted as instantaneous frequency. This can be useful in a number of applications that require the extraction of such parameters.


When upsampling, you don't really stretch the signal in time. You insert new samples between the existing ones, without modifying the times at which those samples were taken. One property of upsampling is that the waveform remains exactly the same before and after the process; you just increase the number of samples. Note that the main difference between oversampling and upsampling is that the former occurs at the time of sampling, and the latter occurs after sampling has already been done. If there is no aliasing, in theory both produce the same result.

Aliasing is only present when there is a signal at the input of your sampler whose frequency is higher than the Nyquist frequency. In your example, if you sample a 70 Hz signal at 280 samples per second, you will not have aliasing. However, in practice you will not always have precise knowledge or control over the signal you're sampling. For example, you may want to sample a signal coming from an antenna. You don't know in advance what signals are going to be picked up by the antenna. In a case like this, one approach is to low-pass filter the signal before sampling, to ensure no aliasing will happen.


If by "oversampling", you mean interpolating more frequency samples between your existing frequency samples, then in this case this would have the effect of zero padding your time domain waveform. Therefore the duration of your signal itself doesn't change, but the duration of time over which the DFT is performed would, with the addition of zero padding to fill out the added time duration (which could be before or after the waveform, changing the phase component of the resulting frequency spectrum accordingly).

To understand zero padding it is helpful to understand the difference between the Discrete Fourier Transform (DFT) and the Discrete Time Fourier Transform (DTFT). The Fast Fourier Transform (FFT) is an efficient algorithm for computing the DFT.

Zero padding a discrete waveform in time results in interpolation of DFT result in frequency. The DFT result is discrete, given by samples of the DTFT. The DTFT in contrast is a continuous function in frequency. By zero padding, we introduce more samples on this same DTFT.

Comparing the formulas for the DFT and DTFT will give you further insight into why this is so:

Here the plots show how the DFT is samples of the DTFT, and by adding zeros to the time-domain waveform prior to taking the DFT (zero-padding), we are approximating the continuous DTFT by filling in more samples (interpolation)-- the more zeros we add the closer we can get to the continuous waveform DTFT. So zero padding and interpolating the frequency spectrum are Fourier Transform pairs.

Further notes:  In the plots I say that the DFT is repeating in time, but in actuality it is only given over N samples. The repetition in time is a mathematical equivalence: similar to the Fourier Series Expansion which is defined over a finite time T, it's frequency components only exist at integer multiples of 1/T (discrete in frequency). If you allowed the time domain waveforms of those components to extend to $\pm \infty$, the base waveform that was defined over t=0 to T would repeat in time. The DFT as given above would be the same as the DFT of the same waveform repeating in time, and seeing that helps provide a lot of insight--- anything that repeats in time MUST be discrete in frequency, anything that is discrete in time MUST repeat in frequency (A/D sampling is an example). The DFT does both. Since the DTFT is performed on a single time domain waveform over all time, it never repeats in time and for this reason has a continuous frequency response.


How to Remove the Patch Artifacts of Neural Network Denoising Process?

One simple way to solve it is using Overlapping Patches.

Let's say you have image which is $ 20 \times 20 $ and you work on patches of the size $ 5 \times 5 $.
As I understand from your description you do 16 times denoising of $ 5 \times 5 $ patches.

What you should do is run the patches mask like in convolution. So each pixels (Ignoring boundaries) will be filtered in 25 different patches.

Now you have 25 different estimates for each pixel.
The trivial solution is to use their mean.
More advanced methods will use some prior on how to create a good weighted combination as output.
One prior could be keeping the composition of the image as smooth as possible.
Other methods could be built around estimation of the noise in each patch.


How to synchronize audio and video streaming for not fixed FPS streaming?

Container formats are usually storing video and sound information in chunks that are supposed to be triggered at the same time. In turn, each chunk might further be compressed according to the specifications of a given CODEC.

For example, consider a video at $FPS=30$ with its sound sampled at $Fs = 48kHz$.

In this example, one frame of the video sequence (one video "sample") should remain on the screen for $\frac{1}{30}$ seconds. But, $\frac{1}{30}$ seconds is equivalent to $\frac{48000}{30}=1600$ audio samples. And the idea here is that you put the frame on the screen and then trigger the sound and by the end of that sound clip you put the next frame up and trigger the next clip.

Or, to be more realistic, you can take the two streams and establish two functions that map audio samples to video samples and the reverse.

For instance, which video frame would be shown to the viewer at sound sample $n_{audio}=987600$?

That would be:

$$n_{frame} = \left \lfloor \frac{n_{audio}}{Fs} \cdot FPS \right \rfloor$$

Where $\left \lfloor \cdot \right \rfloor$ denotes the floor function.

And in this case:

$$n_{frame} = \left \lfloor \frac{987600}{48000} \cdot 30 \right \rfloor = 617$$

If either of your video frame rate or audio sampling frequency are variable, then the respective parts of the above equations become cumulative sums. So, to find where in time are you within either sequence, you have to accumulate the rates that brought you to that time (first 20 frames at 12FPS, next 40 frames at 20 FPS, ... and so on and similarly for the audio). CODECs will have that information in timestamps or even better, the decoder would be able to return "raw" streams and then you are back to the simple equations from above.

But in this case, unless you operate the CODEC yourself (even an uncompressed CODEC, just to specify which sound chunks match which video frames), it will be very difficult to "pipe" two variable streams to a third party software and let it figure out how to put things together.

The simplistic solution here would be to assume that your video runs at the highest FPS and handle other rates by repeating (or dropping if required) frames. The codec at the end of this will assume that the frames and sound it is passed are at a constant rate and "chop them up" accordingly. But I don't know how suitable this might be for your application (?)

Hope this helps.


How to interpret frequency response of moving average system in discrete time signal processing?

In discrete time, a frequency of zero is just a constant sequence, in complete analogy with continuous time. The main difference between discrete time and continuous time is that in discrete time you have a maximum frequency, whereas in continuous time you (theoretically) don't have one. That maximum frequency in discrete time is achieved for a sequence with alternating values, i.e., the sequence $(-1)^n$ has a maximum frequency of $\omega=\pi$. There can be no larger frequency.
If the sequence was obtained from sampling a continuous-time signal with a sampling frequency $f_s$, then $\omega=\pi$ translates to an actual frequency in Hertz via $f=\omega/(2\pi f_s)$.
Returning to the example of the discrete-time lowpass filter, from the frequency response you see that higher frequencies (up to the maximum frequency $\omega=\pi$) are attenuated. That means that a discrete-time sinusoid's amplitude is changed according to the magnitude of the filter's frequency response at the sinusoid's frequency.
So if the input signal equals
$$x[n]=A\sin(n\omega_0)\tag{1}$$
the amplitude of the corresponding output signal is given by $A|H(e^{j\omega_0})|$. If $\omega_0$ happens to be a multiple of $2\pi/5$ then the input is completely suppressed because these frequencies are the zeros of the filter's frequency response.


If X are known 3D points, then the translation is given simply by:

$$T=\frac{\sum_{i=1}^{n}{X'_{i}-RX_{i}}}{n}$$

In case of points being in projected coordinates, then more work have to be done. You can express the 3D points in homogenous coordinates:

$$X=(x, y, z, w)$$

where $w$ is 1 by default.

Then both the translation and rotation can be joined via matrix multiplication:

$$X'=TRX=\begin{bmatrix}1 & 0 & 0 & t_{x} \\ 0 & 1 & 0 & t_{y} \\ 0 & 0 & 1 & t_{z} \\ 0 & 0 & 0 & 1 \end{bmatrix}RX$$

Furthermore, you would be able to add camrea projection matrix in the equation:

$$X'=TRPX$$

Now each point correspondence ($X\leftrightarrow X')$ will give you 4 equations. The minimum number of point correspondences depends on the number of unknowns.

There are usually more than enough points and the problem is overdetermined. The above equation can then be solved in a least-squares fashion.

Another approach is to first perform camera calibration to determine matrix $P$. Then the problem will be simplified greatly.

The best reference I know is (what else) H&Z book. Chapter 7 is all about computing the camera matrix.

Estimation methods are described as well.


Known problems with IIR DF2 when using 64-bit floating point processing?

Depends on how close your poles are to the unit circle. For DF2, the transfer function from the input to the state variables are given by the poles, so there can be substantial gains in there. For typical audio filters, it's not unusual to see 100+ dB of gain and that would definitely be a problem for 32-bit floating point.

These enormous gains also make updating the filters dynamically very tricky since it will typically generate massive transients.

I suggest looking at transposed form 2. It's often as efficient as DF2 but has  much better behaved state variables.


Fit a line pattern on curve with unknown number of points

A useful approach could be to take the numerical derivative of the function. Then you will get some constant line at some $y=a$, where $a$ is non-zero. This would help to eliminate all the other parts and to determine the starting and eding points of the line.

About the fit Matlab has a polyfit command. If the curve of interest is linear than basically this command with parameter $2$ as input will give you the pair $(a,b)$ s.t. the best $ax+b$ fits into the numberical numbers that you provide.


calculating similarity scores between similar and different classes?

From a general perspective, you are trying to take SSIM - which is image difference measure. But it is probably a bad idea if you are looking for measure the intraclass, distance.

You can look at distance measures like Euclidean distance, Mahalonbis distance (if you have co-variance metrics) or others.

As a very simple example : To compute Euclidean distance $L_2$, which is computed between the centroid of the Clusters - as $ Distance = \sqrt{(x_a- x_b)^2 + (y_a- y_b)^2 }$

NOTE: here the dimension are NOT image dimension but are rather the dimension of the feature vector.

  I WILL EXPAND THIS ANSWER WITH YOUR SPECIFIC NEEDS.


Convolution sum vs auto-correlation problem

There are two things to consider. First, looking only at the exam question you can say that the convolution of $h[n]$ with a sequence $y[n]$ is given by

$$(h*y)[k]=\sum_nh[n]y[k-n]$$

Now define $y[n]=x[-n]$:

$$(h*y)[k]=\sum_nh[n]x[n-k]$$

Consequently, we have

$$\sum_nh[n]x[n+k]=(h*y)[-k],\quad y[n]=x[-n]$$

So you need to flip the signal $x[n]$, compute the convolution, and flip the result. This is true, given the sum in your exam question. No need to worry about complex conjugates because we just used the definition of the autocorrelation sum in the question.

The other issue is the definition of the autocorrelation for complex sequences, and how it can be implemented by convolution. The usual definition is

$$R_{xx}[k]=\sum_nx^{*}[n]x[n+k]$$

The convolution of $x^*[n]$ with a sequence $h[n]$ is

$$(x^**h)[k]=\sum_nx^*[n]h[k-n]$$

For $h[n]=x[-n]$ we have

$$(x^**h)[k]=\sum_nx^*[n]x[n-k]=R_{xx}[-k]$$

So for complex signals, you need to compute the convolution sum of $x^*[n]$ with $x[-n]$, and flip the result to arrive at the autocorrelation.

Note that $R_{xx}[k]=R_{xx}^*[-k]$, so instead of flipping the result, you can also use its complex conjugate.


Camera and laser scanner calibration final approach questions

The solution I gave in the other answer assumed that the plane in question was fixed relative to the camera location.  Thus a mapping from the 2D pixel locations to a 2D location on the plane could be made.  Using the same calibration points, a reverse mapping can also be made.  These mappings contain the distortion of the lens, the location of the camera relative to the plane, and the perspective effects.  The mappings can be made more accurate by using higher order equations and more calibration points.

Once you have the 2D coordinates on the plane, call them (x,y), you can convert into real world 3D coordinates using a vector parameterization of the plane:

$$ \vec r =  x \vec a + y \vec b + \vec c  $$

Where $\vec a$ is the 3D real world unit vector in the x direction on the plane, $\vec b$ is the 3D real world unit vector in the y direction on the plane, and $\vec c$ is the 3D real world unit vector of the origin on the plane.

If the situation is more complicated than that, e.g. the laser plane moves, or the camera moves, then a different solution is called for.

I have Python code for a generalized mapping solution.  You can contact me at the email address on my profile page.

Hope this helps.

Ced


what does it mean to have a decorrelated colour space?

edit: to be clear this answer describes why Lab can be described as a decorrelated color space. This does not imply that decorrelation is the main benefit of using Lab (see many answers on why Lab is useful)

If you plot all the RGB colors of a standard RGB image of the natural world you will notice something (see below), the values tend to fall on a diagonal line that isn't one of the axes of the color space.

In fact if you know the red and green values of a color in RGB color space you could make a good guess as to the blue color, it would be the blue that makes the point fall near that clustered line. This is because the values are all very correlated, they are good predictors of each other.

  As I conceive it, CIE Lab is decorrelated means that there's an independence between channels, if for example the value of l* grows bigger, this doesn't necessary imply that the values of a* or b* increase too, but so is the case in rgb I guess. So, what am I missing

So as you can now see, as the values of one RGB dimension increase, the others are also likely to increase.

As it turns out, that clustered diagonal line running through the RGB cube represents the luminance dimension.

In CIE Lab, that luminance dimension we see in the RGB cube as a diagonal running from (0,0,0) to (1,1,1) has been made into one of the axes of the color space. If you know the a* and b* values of a color you still don't know much about the L value. The CIE Lab dimensions are decorrelated from each other because they don't predict each other.


Maximum frequency of any mobile phone speaker

The native sampling frequency of iOS devices is 44100 Hz, while Android devices can operate at either 8, 16, 24, 44.1 or 48 kHz. Support of all of these depends on Android device manufacturer, however most Android support all of them.


Convolution engine for FIR and IIR equal?

Is there any difference between convolution engines for IIR and FIR filters?
  And is there a difference in sound quality using either convolution engine?

A convolution engine performs convolution the discrete version of which is:

$$y[n] = \sum_{k=0}^{|h|-1}x[n-k] \cdot h[k]$$

Where:

$n$ denotes discrete time
$|\cdot|$ denotes length of sequence (so $|h|$ means, the length of the impulse response).
$y$ is the output
$x$ is the input
$h$ is impulse response

Notice here that this formula is giving you $y[n]$. This means that to calculate one time instance for the output you need to take into account $|h|-1$ time instances in the past of the signal ($x[n-k]$ for $k=0$ is the present value).

Both Finite Impulse Response (FIR) and Infinite Impulse Response (IIR) filters can theoretically be expressed via convolution.

But, in the case of IIR filters, $|h|=\infty$. Therefore, to produce a given $y[n]$, you would have to wait forever. What is really done in practice given an IIR filter is to truncate its impulse response when this has decayed sufficiently. The alternative (still in the case of IIR filters) is to implement them in one of the standard forms. In that case, you are implementing a recursive difference equation that produces one $y[n]$ as a combination of a finite set of past values of both its input and output. As a side note, the fact that the filter now has feedback from its output back to its input is what makes it IIR as there are now conditions by which numbers can circulate forever in that loop without necessarily decaying to zero.

Therefore:

When referring to a "...convolution engine" we are talking about FIR filters.
Because of the fact that you would have to truncate the IIR impulse response to express it via convolution, yes the result can be different.

  I could not find much info about the convolution engine within Equalizer APO. I‘m just interested in best audio quality.

Equalizer APO can do both types of filters (look into filters/) plus special implementations for some standard filters. For example, if you look into filters/IIRFilter.cpp and specifically the function IIRFilter::process(), you will see exactly the two iterations that apply coefficients to an input buffer, an output buffer and finally produce an output value. And if you look into filters/ConvolutionFilter.cpp you will notice that it performs convolution via FFT (faster) and libHybridConv (see libHybridConv/libHybridConv.c, hcProcessSingle()).

In conclusion: "Best audio quality..." depends on the workload of your machine. A sufficiently long IIR filter would not sound too different than an FIR one. And at the same time, long filters take more time to calculate. This is further worsened if you want to apply the filter in more than one channels (and I am not talking multiple speakers here, just going stereo means applying one convolution per channel. You should also take into account the comments about the sound quality being subjective here :) ).

However, you might be able to take the results of your room's impulse response, prioritise which "problems" you want to correct and then use a simple graphic equaliser to correct those problems in particular.

You can run a simple graphic EQ, whether in Equalizer APO or any other VST host and then re-route your audio source to be going first through the VST filter and then to the output.

Hope this helps.


What is the impulse response of the system with this difference equation:

Starting from the difference equation
$$y(k)−\frac{5}{12}y(k−1)+\frac{1}{24}y(k−2)=x(k)−\frac{1}{2}x(k−1),$$
you can express the $Z$-transform on both sides ($X(Z)$ and $Y(Z)$ denote the $Z$-transforms of $x$ and $y$) and get:
$$Y(Z)\left(1−\frac{5}{12}Z^{-1}+\frac{1}{24}Z^{-2}\right)=X(Z)\left(1−\frac{1}{2}Z^{-1}\right).$$

Details on this procedure are quite well explained in $Z$ Transform of Difference Equations.

The transfert function $H(Z)$ is given by the ratio of the output by the input:
$$H(Z) = Y(Z)/X(Z) = \frac{1−\frac{1}{2}Z^{-1}}{1−\frac{5}{12}Z^{-1}+\frac{1}{24}Z^{-2}}.$$
As the denominator is (now) a simple polynomial of degree 2 in $Z^{-1}$, you can quite easily expand the rational fraction $H(Z)$ in $Z^{-1}$ to obtain a series $\sum_{k=0}^{+\infty} h(k)Z^{-k}$ by using formal expressions like:
$$\frac{1}{1-X} = 1+X^{-1}+X^{-2}\ldots$$
and
$$\frac{1}{1+aX+bX^2} = c/(X-r_0)+d/(X-r_1)$$
with $r_0$ and $r_1$ the roots of the degree-2 polynomial $1+aX+bX^2$.

With some computations with the numerator, you should be able to retrieve the solution to your expression, with simple formulas for your $h(k)$, providing you with the impulse response you have looking for.

You can find additional material in Partial fraction decomposition. If you prefer a more self-contained, note-like document, I would suggest the following Difference equation solutions from impulse response, which provides an alternative explanation.


sinus generation and it's fft in radians

Why on first plot fft shows that cosine has a frequency 11 radians?

The array index number of the cosine peak is 11.  The frequency is different.

Matlab indexes from 1, so the FFT bin number will be 10 (as we have a 0 bin in FFTs).

The frequency of that bin will be
$$
\frac{2 \pi (11 - 1) }{N } = \frac{2 \pi 10}{ 198} = \frac{10\pi}{99}
$$
assuming your FFT length was 198.

Note that this will be a normalized frequency between 0 and $2\pi$ (or, if you prefer, $-\pi$ and +$\pi$).

  How to shift fft to achieve frequency domain around 0?

Your use of fftshift does this in the second plot.  However, you need to set the $x$ axis by doing something like:

N = length(sinefft);
delta_omega = 2*pi/N; % rad/sec between FFT bins
omega = -pi:delta_omega:(pi - delta_omega)
plot(omega, abs(fftshift(sinefft)));

Minor Correction I previously had this frequency variable as f, it's actually omega.

  How to calculate now what frequency in Hz has this signal?

So the normalized frequency in radian per second is omega.  That ranges from $-\pi$ to $+\pi$ (modulo some end effects).

The normalized frequency in "Hz" is f. This will be omega/2/pi. This ranges from -0.5 to +0.5.

The frequency in terms of the sampling frequency will be f_true = omega / pi / 2 * fs (where fs is the sampling frequency).


KNN method for different values of K on a handwritten digit data set

The k-NN algorithm all by itself is very simple, and can be used on vectors of any dimensionality to classify them. Let me write it down here for the sake of answer completeness:

store all your vectors from the training set
store the class of each vector together with them
when a new, unclassified, vector arrives, compare it to all the vectors from the training set
select the $k$ neighbours that have the smallest distance to the new neighbour
the class assigned to the new vector is the most common class among those $k$ closest vectors

Now, that said, the value of $k$ you choose depends greatly on the size of your training set. E.g. if you have a training set of size $100$, and set $k=100$, the result will always be the same, equal to the most frequent class in your training set. If you set $k=1$, the new vector will be classified as the geometrically closest vector (Euclidean distance) every time.

If you're just supposed to present the results for different $k$ for the sake of discussion and learning, I suggest you start from $k=1$ and test for all the $k$'s sequentially: it will make a nice graph and a nice analysis. Maybe you can skip the even values of $k$ (for $k=2$, if you get two different classes as the closest, you're unable to make a decision).

What's really important is how are you going to shape your descriptor vectors. Each component in your descriptor vector is a feature you describe numerically which should be able to help you make a decision. You also have the curse of dimensionality to worry about: more components your vector has, slower the decision making process is.

You have two approaches here:

Decide carefully on features that are confirmed to be good at the task (do some research), and limit yourself to a certain number.
Work with features you think would be good, but don't worry too much about how good they are. Then, when you have a big number of features, use a method like Principal Component Analysis (PCA) or Fisher's linear discriminant (FLD) to reduce the dimensionality of your descriptor vectors before the k-NN

Lastly, let me give you some ideas about features you might want to use to form your feature vectors. I'm not guaranteeing for any of them, I'm just writing ideas off the top of my head:

Down-sample (shrink) your sample images to $3 \times 3$ or $3 \times 4$ pixel black and white images, and take the 0/1 values of all the 9 or 12 pixels as features. You might want to do a morphological thining or get the skeleton of the pen stroke before you shrink your images.
Instead of taking each pixel from a down-sampled image, you could take only the count of background/foreground pixels in each row and column in the downsampled image, making for a total of 6 or 7 non-binary features as opposed to 9 or 12 binary ones.
Count the number of holes in your (binarized) image. By "hole", I mean a connected component of background pixels not connected to the image border.
Take a simple measure such as describing the elongation of the object, e.g. the horizontal distance between two horizontally furthest pixels divided by the vertical distance between two vertically furthest pixels.
Use the Hu set of invariant image moments (wikipedia, article) It's very old work, but I've used them while making a letter-recognition project during my studies, and for simple uses such as this they can still be helpfull.


Texture mapping on a 3D face from 2D face image

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

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


Solution to Cross-correlation problem of 2 audio signals

Having looked at how xcorr is defined I think I know what is happening.

There are two issues:

The range of xcorr is $2N-1$ (N = length of signal X). The initial values are for negative shifts and the second half is positive shifts. Therefore subtract 120000 to get the correct shift.
xcorr pads Y with zeros to the same length as X. The first 20000 points will be sample Y and the rest of the data zeros. Therefore the value of the maximum correlation will be the start point, not the mid point of the signal (as I originally thought). To prove this think what would happen if Y was the beginning of signal X.

If we now do this on your data we get a start point of (172000-12000)/8000=6.5seconds as expected. The 0.0001 sec comes from how you've defined the first point in time.

However, this is a pretty bad way to do it as any noise in the signal Y will cause your algorithm to fail. I suspect the aim of this task is to manually right code to calculate the correlation function in a similar way to corr does and then find the maximum.


Why is N-point DFT approximated by the sinc function?

When a waveform of arbitrary length is cut into a length that can be fed to an N-point DFT (with N being finite), that is the same as multiplying that original waveform by a rectangular window of that given finite length.  Multiplication in the time domain is the same as convolution in the frequency domain with the Fourier transform of that multiplying function, which happens to be a Sinc function for a rectangular window (or more precisely, for finite length real data, a conjugate mirrored Dirichlet (or periodic Sinc) function, which is pretty close to a regular Sinc, except near DC or Fs/2 or with very short DFTs.  More details here ).

Another way to look at it is that any pure sinusoid of a frequency that is not exactly equal to one of the DFT basis vectors will thus not be orthogonal to any of the DFT basis vectors.  Thus any non-zero magnitude of that sinusoid must end up partially represented in every element of a DFT result (since only orthogonal vector dot products will end up zero).  To interpolate the magnitude of that original non-orthogonal pure sinusoid, one might use the Whittaker–Shannon interpolation formula which happens to use a Sinc kernel to weight the contributions from all the DFT result elements into which energy from the original sinusoid was spread.


Yes, FIR filters are linear.  No, the length of the impulse response does not have anything to do with whether it is linear or not.  For instance, IIR filters are also linear.

Linearity is determined by whether the filter obeys the superposition principle, nothing else.  So, FIR and IIR filters are linear, median filters are not.

EDIT: Technically Dilip is correct.  Finite and infinite impulse response filters can be non-linear if their impulse responses incorporate non-linear operations, such as squaring.  I was going by the "industry" meaning of FIR and IIR, which are linear because they only use multiplication and adding, but in the mathematical sense FIR and IIR filters can be non-linear.


You'd avoid having this problem by detecting whether FM is present before demodulating it.
The reason is simple: FM demodulation throws away a lot of the info that's still in the RF signal – most prominently, the actual received power – and that info is what you'd want to use to squelch.
That can (unless there's other types of transmissions than audio FM going on on the same channels) simply be done via power detection. There's other, elegant, methods working on the radio signal, but many of them depend on the way you demodulate FM.
If you're really stuck with doing it on the audio side: The output of FM-demodulating radio noise only doesn't have the same statistical properties as audio. Sadly, it mostly depends on the type of demodulation you do, again, but what should work relatively OK would be something like estimating the variance of your signal in multiple subbands. If you've got a similar variance all over the bands, that would be unusual for speech. Music, however, might still be caught with that, so be a bit careful.
TL;DR: Do the presence detection on the radio signal, not the demodulated audio.


Comparing with FIR, why does IIR always gives a bias error, when calculating signal DC?

i think this is about a well-known limit cycle problem with IIR filters, particularly implemented using fixed-point arithmetic.  this is something i wrote a long time ago. from which approximately half of this article is derived from.

maybe later i will copy my thing over here and dress up the equations nicely with $\LaTeX$.


In order to use fdesign.comb you have to do something like the following to specify the parameters of the filter

D = fdesign.comb('notch','N,BW',8,.01);

but then you must generate the filter as you mentioned in your question.

H = design(D); % H is a DFILT

you can then view the filter like this

fvtool(H)

One way to apply the filter to your data is to use the filter function

y=filter(H.Numerator,H.Denominator,inArray)

My favorite way to design a filter I've never used is with fdatool because you can play around with all of your parameters. Here is how I would do it:

Type fdatool at the command prompt and a window pops up
Specify Notching in the dropdown that initially says Differentiator
Inder IIR choose comb
Input your parameters
click Design Filter

Once you play around and it looks as you wish, you can either export the coefficients or generate the code that would create this filter.

To export coefficients

go to file and click export
make note of the variable names and click export
now you can use the filter command with y = filter(num,den,data);

To generate code

go to file and click generate MATLAB code
save the file

You see in this file something like the following

Fs = 48000;  % Sampling Frequency

N     = 10;    % Order
BW    = 1200;  % Bandwidth
Apass = 1;     % Bandwidth Attenuation

[b, a] = iircomb(N, BW/(Fs/2), Apass);
Hd     = dfilt.df2(b, a);

the line with iircomb is what generates the filter coefficients, you can now use y=filter(b,a,data);

Hope this helps.


Is there anything called uncorrelated filter?

I do not know if you are referring to a Whitening filter (check this). A Whitening filter can be viewed as a tool converting the correlated sequence into a white sequence. The application of this can be found in Wiener filter and speech compression. For example, in speech compression, the speech signals are highly correlated, which means they are highly redundant. Whitening them by converting them into white sequences (No redundancy) can reduce the complexity to compress the speech.


100 kHz sinewave contaminating a signal sampled at 48 kHz

In your case, since 100kHz is more than twice 48 kHz, it will be folded twice in the sampled version, then there are two pairs of positions: +/-4 KHz and +/- 44Khz.


Should i high pass audio for real time frequency detection?

Yes, it is. If you pickup signals with a microphone you will often see enormous amounts of low frequency signals (air condition, HVAC, traffic, wind, etc.). Even in quiet rooms there tends to be a lot of very low frequency noise.

This noise can dominate the time domain waveform and reducing it can often make your pitch detector more reliable and/or allows for a cheaper or more efficient algorithm.


Calculating an incoherence property from sub-optimal sampling patterns

$\Phi$ is the matrix that represents the way you sample your signal $x$. Actually,  $\Phi$ can be an identity matrix with some rows eliminated, which means you are picking a subset of $x$.

$\Psi$ represents the basis that you choose to expand $x$:
$$x = \Psi c$$
A very simple example is that $x$ is a cosine function (thus containing only one frequency component), then you can take $\Psi$ to be the Discrete Fourier (or discrete cosine) matrix and $c$ is just the coefficient of the specific frequency, and it is pretty sparse (all the coefficients on other frequencies will be zero).

Then the compressed signal $b$ can be represented as:
$$b = \Phi \Psi c$$

Actually you can also view $\Phi \Psi$ as selection of the random rows of $\Psi$, then
one possible option of $\Phi$ and $\Psi$ is (in Matlab):

Psi = fft(eye(n));
p = randperm(n);
Phsi = Psi(p(1:m),:); % choose the first m rows from the permutation

Psi is $\Psi$, and Phsi is $\Phi \Psi$.

For your case, yes you can use singular value decomposition on the image space. But I would prefer to

[u,s,v]=svd(M'*M); % spectrogram
Psi = s(1:m,1:m) * v(1:m,:);
U = u(1:m,:);

which indicates that you select the m most magnitude on spectrum components. Columns in v' are the orthogonal vectors that map the frequency component with the corresponding magnitude to spectrum space. Please note that you don't need u as the basis, since it is always equivalent between u * s * v' * x = b, and s * v' * x = u' * b. As a result, the compressed form is converted to:
$$U'b = \Phi \Psi c$$

Hope it helps. Thanks


i think i found it

X= [x1, x2]
Y=[y1, y2]

c = [[1; 1]  X(:)]\Y(:);
slope_m = c(2)
intercept_b = c(1)

aLine = [slope_m,-1,intercept_b];
points = lineToBorderPoints(aLine,size(negfillHoles))
line(points([1 ,3]),points([2 ,4]),'LineWidth',1.5,'Color','r');
 plot(x1,y1,'*','LineWidth',2,'Color','g');


Your procedure is correct. But, you need to consider two things:

(i) Pole/zero locations.

(ii) Whether the system is causal or non-causal.

Your final answer is OK provided you consider that all the zeros and poles are distinct and the underlying LTI system is causal.

You will get two more solutions if the LTI system is considered non-causal, even if you assume to have distinct zeros and poles.

Hope this helps.


The ffmpeg implementation of the HEVC decoding follows this formula instead:

$ |p_2,_0 - 2p_1,_0+p_0,_0|+|p_2,_3 - 2p_1,_3+p_0,_3|+|q_2,_0 - 2q_1,_0+q_0,_0|+|q_2,_3 - 2q_1,_3+q_0,_3|<\beta \ (1)$

This appears in the hevc_deblock.asm file.

;compare
    pcmpgtw         m15, m13, m14
    movmskps        r13, m15 ;filtering mask 0d0 + 0d3 < beta0 (bit 2 or 3) , 1d0 + 1d3 < beta1 (bit 0 or 1)
    test            r13, r13
    je              .bypassluma

So I conclude that there is a mistake in the paper and the correct formula is that one (since it also makes more sense).

edit: As the comment suggested, I checked in the reference implementation for the algorithm and the interesting part is here in the file TComLoopFilter.cpp:

for (UInt iBlkIdx = 0; iBlkIdx<uiBlocksInPart; iBlkIdx ++)
      {
        Int dp0 = xCalcDP( piTmpSrc+iSrcStep*(iIdx*uiPelsInPart+iBlkIdx*4+0), iOffset);
        Int dq0 = xCalcDQ( piTmpSrc+iSrcStep*(iIdx*uiPelsInPart+iBlkIdx*4+0), iOffset);
        Int dp3 = xCalcDP( piTmpSrc+iSrcStep*(iIdx*uiPelsInPart+iBlkIdx*4+3), iOffset);
        Int dq3 = xCalcDQ( piTmpSrc+iSrcStep*(iIdx*uiPelsInPart+iBlkIdx*4+3), iOffset);
        Int d0 = dp0 + dq0;
        Int d3 = dp3 + dq3;

        Int dp = dp0 + dp3;
        Int dq = dq0 + dq3;
        Int d =  d0 + d3;

//...

        if (d < iBeta) //<-not a >Beta as the paper said
        {
          Bool bFilterP = (dp < iSideThreshold);
          Bool bFilterQ = (dq < iSideThreshold);

          Bool sw =  xUseStrongFiltering( iOffset, 2*d0, iBeta, iTc, piTmpSrc+iSrcStep*(iIdx*uiPelsInPart+iBlkIdx*4+0))
          && xUseStrongFiltering( iOffset, 2*d3, iBeta, iTc, piTmpSrc+iSrcStep*(iIdx*uiPelsInPart+iBlkIdx*4+3));

          for ( Int i = 0; i < DEBLOCK_SMALLEST_BLOCK/2; i++)
          {
            xPelFilterLuma( piTmpSrc+iSrcStep*(iIdx*uiPelsInPart+iBlkIdx*4+i), iOffset, iTc, sw, bPartPNoFilter, bPartQNoFilter, iThrCut, bFilterP, bFilterQ, bitDepthLuma);
          }
        }


What are the eigenvalues of the 8 point DFT matrix?

The eigenvalues belong to the same set of quartic roots of unity verifying $\lambda^4=1$, whatever the order of the DFT.

For more details on their multiplicity, you can read:
Eigenvectors and Functions of the Discrete Fourier Transform, 1982, Dickinson and Steiglitz (online).


Cross-correlation of filtered random processes

As an addition to Dilip's answer I'll show you how to derive that result:

$$\begin{align}R_{y_1,y_2}(\tau)&=E[y_1(t+\tau)y_2(t)]\\&=E\left[\int_{-\infty}^{\infty}x(\alpha)h_1(t+\tau-\alpha)d\alpha\int_{-\infty}^{\infty}x(\beta)h_2(t-\beta)d\beta\right]\\&=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}E[x(\alpha)x(\beta)]h_1(t+\tau-\alpha)h_2(t-\beta)d\alpha d\beta\\&=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}R_x(\alpha-\beta)h_1(t+\tau-\alpha)h_2(t-\beta)d\alpha d\beta\\&\stackrel{\gamma=\alpha-\beta}{=}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}R_x(\gamma)h_1(t+\tau-\alpha)h_2(t-\alpha+\gamma)d\alpha d\gamma\\&\stackrel{\zeta=\alpha-t}{=}\int_{-\infty}^{\infty}\underbrace{\int_{-\infty}^{\infty}R_x(\gamma)h_2(\gamma-\zeta)d\gamma}_{(R_x\star h_2^-)(\zeta)}\; h_1(\tau-\zeta) d\zeta\\&=(R_x\star h_1\star h_2^-)(\tau)\qquad\qquad\qquad (1)\end{align}$$

with $h_2^-(t)=h_2(-t)$.

The cross-spectral density $S_{y_1,y_2}(\omega)$ is the Fourier transform of the cross-correlation function:

$$S_{y_1,y_2}(j\omega)=S_x(j\omega)H_1(j\omega)H_2^*(j\omega)\tag{2}$$

where $S_x(j\omega)$ is the power spectral density of $x(t)$, and $H_1(j\omega)$ and $H_2(j\omega)$ are the frequency responses of the two filters.

Using $(2)$ it is straightforward to define a condition on $H_1(j\omega)$ and $H_2(j\omega)$ such that the cross-spectral density, and, consequently, the cross-correlation function become zero.


Image Gradient - Flipping Effect of the Convolution Matrix (Kernel) for Edge Detection

This is a common question and there are a few other answers on this site which may be helpful. What you have noticed is correct. To be clear, convolution is defined by this flipping, it gives the convolution operation some nice maths properties. Non-flipping convolution is called cross-correlation.

For these kernels, the only thing flipping changes is the sign of the result. You can safely use the inverse (flip) of the kernels to get the form you desire.


Camera Pose Estimation from Vanishing Points

You can use the vanishing points to calibrate the camera.
I recommend reading  Chapter 8.6 of Multiple View Geometry in Computer Vision [ http://www.cambridge.org/9780521540513 ]

Even better, if you have those points on a known sized plane, then you can compute an homography H and the it is possible to compute both intrinsic and extrinsic camera matrices [K],[R|T] respectively.

The size of the object helps determining the scene scale, or more precisely the scale of the column vector t in the extrinsic matrix.
(Zhang's calibration method)[ https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/tr98-71.pdf ]


IIR Filter structure implementation dynamic range

Peter's correct.  at least if fixed-point arithmetic is used.  the DF2 has poles before zeros.  that means the signal is getting boosted by the poles before the zeros (which are often sitting very close to the poles) beat the gain back down to reasonable levels.  so if the signal overflows, saturates, and distorts due to this gain, the attenuation offered by the zeros does not reverse that distortion error.  it's too late.

now the DF1 implements the zeros before the poles, and if you did not have an accumulator that has a double word width, then the DF1 suffers a sorta opposite flaw.  instead of potentially overflowing (and saturating) due to gain from the poles which is what the DF2 does, the DF1 will have attenuation first (due to the zeros) followed by the compensating gain from the poles.  if the accumulator is no wider than the signal word going in, there is roundoff error due to quantizing the attenuated signal.  then when the poles act on the attenuated and quantized signal, that quantization error is boosted by the poles.

BUT if your DSP (or whatever processor you use) has an accumulator that is twice the width of the two multiplicands, there is no roundoff error for the poles to boost.  so the DF1 avoids the overflow problems of the DF2, but if and only if there is a double-wide accumulator, does not suffer from boosting the roundoff error.


Difference between these frequency response curves

The two frequency response plots are the result of two different models. In the first one, the model is a mass-spring-damper system as shown here, and the plot shows the ratio of the amplitudes of mass displacement and driving force.

The second plot shows the transmissibility of a single DOF mechanical system as shown below:

(from Fundamentals of Mechanical Vibrations by S.G.Kelly).

Note that here, unlike in the other model, the support is not fixed. The transmissibility shown in your second plot is the ratio of the amplitudes of the mass displacement $x(t)$ and the (forced) displacement of the support $y(t)$.

In formulas, your first plot shows the function

$$M(\omega)=\frac{1}{\sqrt{\left(1-\left(\frac{\omega}{\omega_0}\right)^2\right)^2+\left(2\zeta\frac{\omega}{\omega_0}\right)^2}}\tag{1}$$

where $\omega_0$ is the system's natural frequency, and $\zeta$ is the damping ratio. In electrical engineering terms, this is the magnitude of the frequency response of a second-order low pass filter.

The second plot shows the transmissibility given by Eq. (3.73) in Kelly's book:

$$T(\omega)=\sqrt{\frac{1+\left(2\zeta\frac{\omega}{\omega_0}\right)^2}{{\left(1-\left(\frac{\omega}{\omega_0}\right)^2\right)^2+\left(2\zeta\frac{\omega}{\omega_0}\right)^2}}}\tag{2}$$

Note that $(1)$ and $(2)$ have the same denominator but different numerators.


Perhaps transferable to your problem, one method from Olyha B., Rutledge J., "TrackPoint System Version 4.0 Engineering Specification", IBM, 1999:

  2.2.3 DRIFT CORRECTION

  Due to the significant temperature sensitivity of the force sensors, it is necessary to recalibrate the zero force
  origin of each axis on a periodic basis. In order to properly
  recalibrate, it is necessary to determine when a ‘hands off’ condition
  exists (i.e. when no force is being applied). This is done by
  monitoring the X and Y axis force samples, and when there is minimal
  change in the readings (no more than ±1 count) for a specified length
  of time (3 seconds by default) it is assumed that the hands off
  condition exists. At this time the X, Y, and Z axis origins are reset
  to the smoothed current value of the X, Y, and Z axis force samples
  respectively. It should be noted that noise introduced by any source
  other than applied sensor forces will severely impact the
  effectiveness of this technique.

In my experience as a long-time ThinkPad user, this technique cannot always differentiate between the user trying to move the mouse cursor at constant speed and real sensor drift. Then again this is not usually dangerous and it is easy to manually trigger recalibration by hands off. Perhaps there should be a more advanced analysis algorithm running that can identify velocity variation in a manual attempt to move the mouse cursor at constant speed.


Detection of sharp peaks in an audio signal

"Sharp peaks" in the time domain translates to high frequency content. Your simple algorithm of min and max values would be improved by first passing your data through a high pass filter - otherwise slow (low frequency) rumbles could degrade your detection methods.

To take your solution to the next level, you should use an "onset detection" algorithm. These are well studied. A recommended paper to start is:

A Tutorial on Onset Detection in Music Signals - Juan Pablo Bello et al, September 2005.


Revert anayltic signal(length(N/2)) to real signal (length(N))

First answer deleted since I completely misunderstood the question.

The problem here is that apparently you are missing one frequency point. From the question it's not clear what is missing. You say you have the points from 0 to N/2. That's actually N/2+1 (if 0 and N/2 are included) and all you need. If that's not the case, there are gew possibilities

Go to the person who provided the data and ask them to fix it
X(0) and X(N/2) are actually real valued. Sometimes people pack the X(N/2) value into the imaginary part of X(0) to save space and create a symmetrical representation. Check your X(0) and figure out whether it's imaginary part is actually 0 or not
There is a way to represent a N-point real sequence with exactly N/2 complex coefficients that involves a  half sample rotation in the frequency domain. That's pretty esoteric and unlikely in this case (but works quite well).
If the high frequency energy in your original signal is sufficiently low, you can set the Nyquist coefficients simply to zero. According to the sampling theorem, that should be zero anyway. In practice this will add some noise to the reconstructed signal but if it's small enough that may be tolerable.