How do I normalize accelerometer signals values to a relative scale?

Assumptions

There is a no function in (all) target operating systems that gives you the maximum/minimum values of the accelerometer
You do not want to spend a lot of money and time to build a database of these values on your own

Suggested solutions

I assume there is no way of getting around a kind of calibration procedure. Here are two thoughts:

#1

The simplest one that I can think of is that you let people drop their phone from a certain height, let's say 1 meter. Have them to do this:

Hit a calibrate switch. This starts a count down that starts some kind of user feedback (e.g. a certain color on the screen, a tone being played) for 5 seconds.
Within these 5 seconds the phone is to be held still, approximately 1 meter above a soft ground that is safe to let a phone fall onto (e.g. blanket, bed, ...)
When the 5 seconds are over and the feedback stops (or changes, e.g. different color on the screen, or the tone stops playing), the phone is to be dropped on the soft ground.

Since you know approximately the height of the phone, you can select a part of the accelerometer data where the phone certainly is in free fall, and is therefor free of vibrations from the users hand or from the phone touching the soft surface. All you probably need is to get 10 centimeters of plain, free fall. Rotations are not a problem if you intend to use $\bar{a} = \sqrt{a_x^2 + a_y^2 + a_z^2}$. If needed, you could also let them repeat the process to gain a higher SNR from averaging. For your application, it is certainly sufficient to assume that the measured signal is due to the earth's surface acceleration of $g \approx 9.81\,\frac{m}{s^2}$. If you need slightly higher precision, you could use the geolocation of the device and an theoretical gravity approximation model.

#2

If the app is used in certain contexts and only those values are important, you could simply introduce a learning phase. Let the user wear the device/use it in the intended way, but only for learning. You could:

let the user rate the acceleration during the training phase and use this as a ground truth for your model
just collect the data and let the app figure out the levels by itself, assuming that after let's say... 13 training sessions all intensities were measured

Further ideas

In both cases, you could ask the user to share this data with the other users of your app: After a successful training session, ask them if the would share the phone model/accelerometer model (that you probably can retrieve from the OS) and the calibration values. This way, you could have the app look for calibration values online. If they are available, they could be downloaded and used without the inconvenience for the user of doing the calibration. Additionally, with such a database, you could process and thereby refine the calibration data.


Best reference for object size estimation in an image

The usual way to do this is to form a calibration plate.  One way to form a calibration plate is to place known-shaped "fiducials" in a known-sized array on the surface where you want to do the measurement.

Then, knowing the measurements of the calibration plate, you can count pixels and do a map from pixels in each direction (possibly in each area of the image) to real-world measurements.

You also need to know that your camera will always be in a fixed location relative to the measurement surface, and that the lighting is good enough for you to pick up the external shape of the object.

An example of one possible calibration plate is below.


Whenever you have doubts regarding the properties of the convolution operator, you should resort to its definition.
Let $x_s[n] = x[n+2]$ and $h_s[n] = h[n+2]$. Then:
$$y_4[n] = x_s[n] * h_s[n] = \sum_{k=-\infty}^\infty x_s[k]h_s[n-k]$$
If we replace with the definitions of our discrete functions:
$$y_4[n] = \sum_{k=-\infty}^\infty x[k+2]h[n-k+2]$$
See what happens if we make the change of variables $m = k+2$:
$$y_4[n] = \sum_{m=-\infty}^\infty x[m]h[n-(m-2)+2]=\sum_{m=-\infty}^\infty x[m]h[n+4-m]$$
Does that final expression sound familiar? You were right indeed in your OP, as that is exactly what you thought:
$$y_4[n] = \sum_{m=-\infty}^\infty x[m]h[n+4-m] = y_1[n+4]$$


You are not really looking for a clustering algorithm in the common meaning of the word, but more of a Graph Theory concept.

  So we have a matrix like

That is not a matrix. That is an adjacency list that states nodes adjacent to others.

The "matrix" you might be after is the adjacency matrix. For a simple graph (a graph that does not contain nodes with self-connections), $G(V, E)$ where $V$ is the set of vectors and $E$ is the set of edges, the adjacency matrix $A$ tells you if there is a connection between two vertices in $V$. For example, $a_{2,3}$ would tell you if there is a connection between vertices $v_2, v_3 \in V$.

Applied to your problem, the adjacency matrix looks like this:

$$V = \{ 1,2,3,4,5\}$$
$$ A = \begin{bmatrix}0 1 1 0 0 \\ 1 0 0 0 0 \\ 0 0 0 0 1 \\ 0 0 0 0 0 \\ 0 0 1 0 0\end{bmatrix} $$

A single graph can be composed of a number of components. In your case, node 4 is completely unconnected and it could be considered as a component.

If your graph was undirected, the task of finding connected components would be very simple. But, your graph is directed. This is denoted by the lack of symmetry of $A$ around the main diagonal ($A_{1,3} \ne A_{3,1}$). In this case, you would be looking for strongly connected component algorithms.

A Python module you can use to compute both is Networkx. With functions for connected components and strongly connected components here and here respectively.

This concept of clustering is very interesting but (at first view) quite different than the traditional clustering of data points in some space based on "similarity". See for example this type of clustering and you will see why you need a concept of "distance". In your type of problem, this would be equivalent to having the nodes also report HOW SIMILAR they are to each other. Then, in addition to direction, you would also have a weight (similarity) associated with each edge. In that case, things get a little bit more complicated.

Hope this helps.


What Are Linear and Circular Convolution?

Linear convolution is the basic operation to calculate the output for any linear time invariant system given its input and its impulse response.

Circular convolution is the same thing but considering that the support of the signal is periodic (as in a circle, hence the name).

Most often it is considered because it is a mathematical consequence of the discrete Fourier transform (or discrete Fourier series to be precise):

One of the most efficient ways to implement convolution is by doing multiplication in the frequency.
Sampling in the frequency requires periodicity in the time domain.
However, due to the mathematical properties of the FFT this results in circular convolution.

The method needs to be properly modified so that linear convolution can be done (e.g. overlap-add method).


z-space replaces multiplication with convolution and frequency eigenvectors

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

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

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

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

EDIT: Lets use a specific example

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

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

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


FFT of High Frequency and Low Frequency?

The frequency span of an FFT is determined by the number of samples fed into it and the sample rate of those samples.  Basically, the maximum frequency analyzed by the FFT will be half the sample rate.  The lowest frequency will be DC, and the remaining output "buckets" will be evenly spaced in frequency up to the (last) max frequency bucket.  The number of output "buckets", in the standard implementation, is half the number of input samples.

(Actually, the total number of outputs is equal to the number of input samples, but half the outputs are the imaginary half of complex numbers ... and at about this point my head starts to hurt.)


Trouble with understanding moving average filter as an oversampling filter

A moving average over 6 samples will not add 6 bits of precision to the ENOB (equivalent number of bits) as stated in the question. In general, a moving average over $N$ samples will increase resolution by $log2(N)/2$, and for 6 samples this would be approximately 1.3 bits.

This approximation holds to the extent that the quantization noise on each sample is represented by a uniformly distributed white noise process (which holds well when there are a sufficient number of quantization levels and the sampling clock is not correlated to the signal being sampled). From statistics, when you add N independent random variables, the standard deviation increases at $\sqrt{N}$. Thus in the average where we also divide that sum by $N$ the standard deviation is scaled by $\sqrt{N}/N = 1/\sqrt{N}$.

This results in a 3 dB reduction in quantization noise for every doubling of samples in the average (from $20Log_{10}(1/N)$ since dB is $20Log_{10}(\sigma)$ where $\sigma$ is the standard deviation. As already detailed in the post linked below, each bit of precision represents a 6 dB reduction in quantization noise:

What are advantages of having higher sampling rate of a signal?

A moving average over 4 samples would add 1 bit of precision, 16 samples would add 2 bits, etc.

To support the increased number of bits, the result from the addition of the samples must have that many number of bits which is also apparent to eliminate overflow conditions. Filtering a digital signal (which a moving average is) grows the signal! This is why we must use extended precision accumulators when filtering. We can scale the accumulated signal back after summing (for the average) while scaling the signal first prior to the summation, or equivalently scaling the coefficients, would result in further noise contribution.


Sampling a continuous signal $x(t)$ can be modeled by multiplication with a delta pulse train:

$$x_s(t)=x(t)\cdot\sum_{n=-\infty}^{\infty}\delta(t-nT)=
\sum_{n=-\infty}^{\infty}x(nT)\delta(t-nT)$$

where $T$ is the sampling interval: $T=1/f_s$ ($f_s$ is the sampling frequency). Multiplication in the time-domain corresponds to convolution in the frequency domain. So the spectrum of the sampled signal can be computed by convolving the transforms of $x(t)$ and of the pulse train:

$$x(t)\Longleftrightarrow X(f)\\
\sum_{n=-\infty}^{\infty}\delta(t-nT)\Longleftrightarrow
\frac{1}{T}\sum_{k=-\infty}^{\infty}\delta(f-kf_s)
$$

Convolution of these two Fourier transforms gives

$$x_s(t)\Longleftrightarrow X_s(f)=
\frac{1}{T}\sum_{k=-\infty}^{\infty}X(f-kf_s)\tag{1}$$

Equation (1) shows that the spectrum of a sampled signal is composed of shifted versions of the original spectrum. The shift corresponds to multiples of the sampling frequency $f_s$. Consequently, if the spectrum $X(f)$ is exactly zero for frequencies greater than $f_s/2$, then the shifted versions do not overlap, and the original spectrum can be reconstructed by lowpass filtering the spectrum of the sampled signal. If $X(f)$ is not band-limited to $f_s/2$, i.e. it has frequency components greater than $f_s/2$, the shifted spectra in Equation (1) overlap, and the original signal can in general not be reconstructed from the sampled signal. This is what the sampling theorem states.

EDIT: For non-periodic signals - i.e. if $X(f)$ does not contain Dirac impulses - it is sufficient for reconstruction that $X(f)$ does not contain frequency components above $f_s/2$. In this case it is (theoretically) possible to reconstruct $x(t)$ from its samples $x(nT)$ by applying an ideal lowpass filter with cut-off frequency $f_s/2$. The impulse response of such an ideal lowpass filter is given by

$$h_{LP}(t)=\frac{\sin (\pi t/T)}{\pi t/T}$$

Filtering $x_s(t)$ with this impulse response yields the original signal

$$x(t)=\sum_{n=-\infty}^{\infty}x(nT)\frac{\sin [\pi (t-nT)/T]}{\pi (t-nT)/T}\tag{2}$$

Equation (2) is usually referred to as the Whittaker-Shannon interpolation formula.

Normally one does not need to consider what is actually happening at the frequency $f_s/2$. The only case where this is necessary is the presence of a pure sinusoid with frequency $f_s/2$. From Equation (1) we can see that in this case the shifted spectra (for $k=\pm 1$) will also produce spectral lines at $\pm f_s/2$, which in general will make impossible the reconstruction of the sinusoid at $f_s/2$. Apart from this somewhat artificial case, it suffices to require that the signal does not contain any frequency components greater than $f_s/2$.

The above considerations are of course of a purely theoretical nature, because there exist neither ideally band-limited signals nor ideal lowpass filters. Consequently, in practice we will choose the sampling frequency such that it is well above twice the maximum signal frequency of interest. The anti-aliasing filter must pass frequencies up to this maximum frequency, then it has a finite transition bandwidth, and then we can choose a frequency in the stopband of the filter where the stopband attenuation guarantees that aliasing errors are sufficiently small for the given application. This frequency can finally be considered the "bandwidth" $B$, and the sampling frequency can now be chosen as $f_s\ge 2B$.


inverse FFT when only a few time points are needed

What you're looking for is called a pruned DFT. In principle, it is possible to calculate a subset of outputs from a DFT using fewer mathematical operations. In practice, however, existing highly-optimized FFT implementations like FFTW are designed for full-output transforms. You'll find in many cases, unless you're only concerned with a very small proportion of the transform outputs, it is fastest to just calculate the entire thing and throw away what you don't want.

As it turns out, FFTW has a page in their documentation that talks about this specifically.

If you want to put in a lot of work, it is possible that you could obtain some performance gain by rolling your own FFT implementation for your specific case. FFTW is a general-purpose, any-size-transform library that works well in a variety of conditions, but it's possible that you could beat it if you tune something specifically for your application. Before expending this effort, I would think long and hard about whether it is truly a limiting factor in your system's performance.


Shift property of the Short-term Fourier transform (STFT)

Of course the STFT changes when the signal is shifted. According to the definition we have

$$\text{STFT}\{x(t)\}=X(\tau,\omega)=\int_{-\infty}^{\infty}x(t)w(t-\tau)e^{-j\omega t}dt\tag{1}$$

If we define a shifted signal $y(t)=x(t-T)$, the corresponding STFT becomes

$$\begin{align}\text{STFT}\{y(t)\}=Y(\tau,\omega)&=\int_{-\infty}^{\infty}x(t-T)w(t-\tau)e^{-j\omega t}dt\\&=\int_{-\infty}^{\infty}x(u)w(t-(\tau-T))e^{-j\omega u}e^{-j\omega T}du\\&=e^{-j\omega T}X(\tau - T,\omega)\tag{2}\end{align}$$

From $(2)$ you see that time shifting results in a phase factor and a shift in the STFT. The magnitudes of the STFTs are related by just a time shift, which is of course identical to the time shift of the signal:

$$|\text{STFT}\{x(t-T)\}|=|Y(\tau,\omega)|=|X(\tau-T,\omega)|\tag{3}$$


When lowpass filtering with a Gaussian kernel, how can $\sigma$ be related to cycles/pixel?

The frequency response of a Gaussian is also a Gaussian, see here. The cutoff frequency $\sigma_f$ is defined as the standard deviation of the frequency response. This is given by $$ \sigma_f = \frac{1}{2\pi\sigma} $$ Given that $N$ is the wavelength of the frequency components you wish to remove, you can calculate the corresponding frequency and work out how much you want to reduce it by, work out the corresponding cut-off frequency and filter standard deviation.

A more back of the envelope approach is to say that Gaussians of standard deviation, $\sigma$, remove or strongly reduce components with lateral size $<\sigma$.


DFT zero-padding of signals starting before n=0

Use the second one Tony... It yields the correct implied phase relationship.


MATLAB's $\tt cpsd$ and $\tt pwelch$ - different results for cross spectral power density

Whilst I have not worked out how to calculate CPSD from power spectra, I have confirmed the following two methods are equivalent (using amplitude spectra, rather than power spectra or cross-correlation):

Fs = 1000;
t = 0:1/Fs:1-1/Fs;
x = cos(2*pi*100*t) + randn(size(t));
y=sin(2*pi*160*t)+cos(2*pi*150*t)+randn(size(t));
N = length(x);
nfft=500;
lfft=2*nfft;
overlap=0;
win=rectwin(N);

xdft1 = fft(x,2*nfft);  <outputs two-sided complex amplitude spectra
ydft1 = fft(y,2*nfft);
Sxy1a=xdft1.*conj(ydft1);
Sxy1b=Sxy1a./N^2;
Sxy1c=Sxy1b(1:lfft/2+1);
Sxy1c(2:end-1)=2*Sxy1c(2:end-1);
Sxy1=abs(Sxy1c)';

%CPSD Method 2
Sxy2a=cpsd(x,y,win,overlap,2*nfft,Fs);   < uses the cpsd function with Welch’s overlapped method
Sxy2=abs(Sxy2a);


Digital Signal Quantized Processor Complexity

First, if there's a more formal way to do this than what I'm going to show, I don't know it.

Call out the MAC operations along with data path widths.  I.e., "the FFT algorithm used takes $4 n \log n$ operations with 24-bit wide data paths"

A corollary to this would be that sometimes you may have different precision in your coefficients and your actual data.  So you may call out "24 bit signal path, 16 bit coefficients".

Call out the LUTs used (in an FPGA).  This isn't direct, and has a lot to do with the cleverness of the implementation, but within a specific FPGA technology it has meaning.

Note that these days, unless you're doing full-custom ASIC design, data paths come in chunks - processors give you multiples of their native word width, and FPGAs come with "DSP blocks" that come in fixed widths.  I'm only guessing, but I'd bet money that if you go to design an ASIC then at the first level of design complexity there will be signal processing libraries with their own fixed-width data paths, and that libraries with variable-width data paths would take more engineering resources to use.


Least Mean Squares (LMS) Filter Weight Update

That really depends on context, but generally adaptive implies that the calculations are done on-line / on the fly.

In some applications, the filter is updated for a while, then the adaptation is turned off and the last lot of weights are used.


Rollercoaster: determining subject position from inertial data?

The ride has temporal variance (ex. sometimes getting from A to B takes 300ms longer)

The carts are driven by a motor on certain parts of the ride (e.g. up a slope). The acceleration should be fairly stable during those parts. Any temporal variation, if not attributed to different acquisition timings, is probably because of aerodynamic factors. For a given ride, that could be temperature of the air (which affects density), number of passengers in the carts (body parts disturbing the flow of air) and direction of wind. With an average speed of 82 mph (~132kmh), the roller coaster advances $\frac{132000m}{3600000ms} \approx 0.03666m/ms$. Therefore, a 300ms variation accounts for a "position error" of $300ms \times 0.03666m/ms = 11m$. Depending on the actual length of the ride, this might be negligible.
In addition to this, you have some jitter on your sampling frequency because the phone's OS doesn't guarantee that it will sample at an absolutely strict frequency.

Phones may have different initial orientations

It doesn't matter, if you subtract the first measurement from subsequent measurements then you have a relative indication rather than an absolute. For example, assuming that the first few measurements of your bearing were [270, 265, 260, 255, . . .] then the series [0, -5, -10, -15, . . .] indicates relative turning.

Phones stay in the subjects' pockets, but jostle a bit

Similarly, compared to the G forces that characterise the turning points in the ride, jostling of a phone in someone's pocket might be negligible. Having your phone tucked snugly in a front pocket is one thing, securing it in a zip pocket of a jacket that might flail about during the ride is another.

Learn the "mean" rollercoaster ride

as a time series
as a sequence of significant rollercoaster events

The "mean" roller coaster ride is the average of all available recordings for a given ride, per sensor. There is nothing to "learn" there, the more recordings you have, the closer your average approximates the "average" experience on the ride because you get a number of different accelleration / rotation / bearing samples for each time instance on the ride.
A sequence of "significant roller coaster events" is easy to extract once you have a clear definition of what a "significant roller coaster event" is.
For example, a "significant roller coaster event" might be a sharp drop, or a tight turn. Both of these are characterised by rapid rates of change in the total accelleration of a body on the cart. Therefore, you can derive a total accelleration from your 3D accelerometer as $A_{total}(n) = \sqrt{x(n)^2 + y(n)^2 +z(n)^2}$ and then apply some "event theshold $T_e$ on your $A_{total}(n)$. The indicator function that results from this (some $Z_n = aT_n \ge T_e$) points to some "significant roller coaster event". In the above, $x,y,z$ are the individual components of the accelerometer's measurements. Obviously, it is also possible to combine more than one sensors to extract a "significant roller coaster event".

On a new ride, in real-time, robustly determine where the subject is along the path so we can align music to their experience \m/
Bonus: 3D model the rollercoaster path

You will not be able to extract an accurate 3D model of the roller coaster path from accelerations because it is impossible to do accurate integration of the phone's sensors to a proper 6 DOF inertial nav system. If you want an accurate model of the roller coaster path, simply record GPS coordinates and you have the shape of the ride in 3D in zero time and with minimal data processing.
Tracking where a user is on a given ride is trivial once you have that "average track" signal generated. With that available, tracking a user on the ride is translated to "given a user's "past history" of a few points, where are they likely to be along this specific one dimensional curve?".
So, given that a "random" phone's signals are expected to vary in approximately a given way (the ride's shape), all you have to do is check where they are. You can do this by establishing a Kalman Filter on "position" which monitors current total accelleration and rate of turn. The Kalman will constantly be making a guess of where it should be along the curve and also be correcting this guess with the actual signals it receives from the sensors. You can achieve the same with a particle filter too, arguably with less processing.
Obviously, you can then set markers on the one dimensional curve to trigger your music replay.
Hope this helps.


Distances In a Single Image With Some Real References

Well, you can have an approximation. Due to the forms in the image are non regular, completely plain, its hard if not impossible to know the camera distortion. But if you know the size of the wheels and also you have the measures of the stones you can approximate the distance of the cart to the stones by using this equation:

$$
D_w = \frac{S_w * f}{P_w},
$$

where: $D_w$ is the distance from the camera to the wheel, $S_w$ is the real height of the wheel, $f$ is the distance of the focal to the film when the film is in the camera (You need to search for this parameter to use this approach, also you need to know the proportion of the real impressed picture with the film) and $P_w$ is the height of the wheel in the film. The same thing for the stone:

$$
D_r = \frac{S_r * f}{P_r},
$$

where: $D_r$ is the distance from the camera to the stone, $S_r$ is the real height of the stone and $P_r$ is the height of the stone in the film. Finally the distance $D_t$ from the wheel to the stone is:

$$D_t= D_r -D_w$$

I recomend you to use the objects that are marked in red in the next image, because are the nearest to the epipolar line.

Hope this could help, Best.


Channel Shortening for underwater acoustic channel communication

When I see "channel shortening," I think of DFE coefficient calculations, regardless of the application. This is used in cases from simple DFE equalization to channel shortening for xDSL with DMT to lessen CP requirements or for RSSE in single-carrier systems. The point is not to actually do a DFE; rather, the point is to use the coefficient calculation and paradigm.

The classical DFE architecture is to have two filters: a feedforward (FF) filter and feedback (FB) filter with a slicer after the FB subtraction and decision-directed feedback passed through the FB filter. The intent is to use a FF filter to compress the residual channel response into a smaller number of taps than you started with. For DFE applications, the intent is also to transform a possibly non-mininum phase response to minimum phase to reduce error propagation. The desired structure is to have a FF filter that is (close to) all-pass and simply transforming phase, while the FB taps of the DFE represent the residual channel. That said, the concept of compressing the channel into a shortened response can generalize to the case where the FF filter is not all-pass. In this case, it looks more like a FF equalizer for a partial-response system.

I suggest investigating DFE coefficient calculation and motivations to see approaches for channel shortening. It may also be instructive to look at partial-response literature. In the case of PR, you may see descriptions of equalization to a predetermined response, but it can also be used more generally for compression of a response into a reduced number of taps.


Role of Power Metric and LPF in AGC(Automatic Gain Controller)

Either power or rms value is typically used for AGC as a peak value would be much more variable given it is an instantaneous measurement. The low pass filter typically refers to the Loop Filter (and if so, would be better described as such as it more typically done as an integrator) which is a core part of the control loop design. A band pass filter prior to the power detector is also important since any signals within the bandwidth of the filter/detector would control the AGC level.

Below are example block diagrams depicting typical AGC implementations in the analog and digital domain.  The log and alog interfaces for the digital implementation provide consistent loop bandwidth over the operational range of the AGC due to using multipliers for the gain adjust.

Further AGC considerations for digital implementation are covered in this post: How to adjust receiver gains to avoid saturation and quantization noise to optimise post digital processing?


How to extract a smooth contour from a set of points in 3D?

From your description it sounds like surface reconstruction methods such as Poisson Reconstruction would fit your description. These algorithms already have the necessary smoother built-in, mainly in the form of a regularizer. To have slightly more accuracy, in general I prefer the following algorithm:

  SSD: Smooth Signed Distance Surface Reconstruction, F. Calakli and G. Taubin, Computer Graphics Forum, Vol. 30, No. 7, 2011, http://mesh.brown.edu/ssd/

There is a nice Adaptive Solvers repository containing many popular methods under a single roof:
https://github.com/mkazhdan/PoissonRecon

I strongly recommend to take a look at this. Here is how you would like to use it:

Export your point set and the surface normals to a 3D file : sample.ply
Run (make sure to play with the depth value):

SSDRecon --in sample.ply --out sample_reconstruction.ply --depth 10 --density

This will generate a reconstruction but if you have an open surface the reconstruction will contain undesired extrapollations. So move to (3). If the surface is perfectly closed, then the output will already be good in that stage.

If you have an open surface then run:

SurfaceTrimmer --in sample_reconstruction.ply --out sample_reconstruction_trimmed.ply --trim 7

Again, make sure to play with the trimming threshold. This will give you very good reconstructions such as (taken from here):

Note that for all these algorithms to work, surface normals are essential. Luckily, they can be computed from the data using local plane fits. Moreover, one also needs to ensure the correct signs of the normal directions to properly orient the surface. In the case that a single camera captures the data, the canonical choice is to flip all the surface normals so that they make an acute angle with the viewing direction.


Rectangular or trapezoidal signal filtering. Beginner

Depending on your noise model, a median filter might be quite good or not. Have a look at this example:

import scipy.signal

# from http://stackoverflow.com/questions/23289976/how-to-find-zero-crossings-with-hysteresis
def hyst(x, th_lo, th_hi, initial = False):
    hi = x >= th_hi
    lo_or_hi = (x <= th_lo) | hi
    ind = np.nonzero(lo_or_hi)[0]
    if not ind.size: # prevent index error if ind is empty
        return np.zeros_like(x, dtype=bool) | initial
    cnt = np.cumsum(lo_or_hi) # from 0 to len(x)
    return np.where(cnt, hi[ind[cnt-1]], initial)

Fs = 300
T = 9
a = np.zeros(T*Fs)
# setup the state changes
a[0*Fs] = 1
a[0.5*Fs] = -2
a[2*Fs] = 2
a[3*Fs] = -2
a[4.5*Fs] = 2
a[7*Fs] = -2
a[8*Fs] = 2
a[8.5*Fs] = -2
# create the rect-waveform
data = np.cumsum(a)

# add some noise to the signal
rx = data + 2*np.random.randn(len(data))

plt.subplot(121)
plt.plot(rx)
plt.plot(data, 'r', lw=2)

plt.subplot(122)
# plot the filtered signal (if you have access to the purple signal in your plot, just use this signal)
filtered = scipy.signal.medfilt(rx, Fs//2+1)
plt.plot(filtered, lw=3)

# Perform thresholding with Hysteresis
detected = 2*hyst(filtered, -0.3, 0.3)-1
plt.plot(detected)

# From the thresholds, get the state changes via the derivative
diffed = np.diff(detected)
print ("estimated times of state switch:", ["%.2f" % (x/Fs) for x in np.nonzero(diffed)[0]])
print ("actual times of state switch:", ["%.2f" % (x/Fs) for x in np.nonzero(a)[0]])

estimated times of state switch: ['0.18', '0.60', '2.05', '3.10', '4.56', '7.07', '8.07', '8.59']
actual times of state switch: ['0.00', '0.50', '2.00', '3.00', '4.50', '7.00', '8.00', '8.50']

It creates an artificial signal, adds some noise, does a recovery of the original signal via median filtering. Eventually, it detects the state changes in the signal.


Relation between Kalman filter and Sequential linear MMSE estimation

The system model for the Kalman filter is:

$
x_k = F_k x_{k-1} + B_k u_{k-1} + w_k \\
z_k = H_k x_k + v_k
$

The system model for sequential linear MMSE estimation might be:

$
x_k = x_{k-1} \\
z_k = H_k x_k + v_k
$

Those are equivalent if $F_k$ is the identity, $B_k=0$ and $w_k=0$ and I believe the estimators turn out to be the same in that case. If $F_k$ is not the identity, or if either $B_k$ or $w_k$ are non-zero, then I see no immediate way to match sequential linear MMSE estimation to the Kalman filter system model.

On the other hand, let's consider the system equations for $k=1$, $k=2$ jointly:

$
z_1 = H_1 F_1 x_0 + F_1 B_1 u_0 + F_1 w_0 + v_1 \\
z_2 = H_2 F_2 (F_1 x_0 + B_1 u_0 + w_0) + F_2 B_2 u_1 + F_2 w_1 + v_2
$

This can be stacked into one equation:

$
z = A x_0 + \tilde{v} \,,
$

where $z$ is a column vector stacking $z_1$ and $z_2$. $A$ contains $H_1, F_1, H_2, F_2$, and $\tilde{v}$ contains $H_2, F_2, B_1, u_0, u_1, w_0, w_1, v_1, v_2$. Both the Kalman filter and the linear MMSE estimator can be applied to this linear model, and I would expect them to produce the same result. The Wikipedia article assumes that the error covariance for $\tilde{v}$ is zero, but we can trivially fix that by subtracting the error mean (which arises from the known inputs $u_0, u_1$).

So if we model the system with joint equations, rather than sequential, iterative equations, then linear MMSE estimation applies directly to the general model, and Kalman filtering gives no improvement for this more general problem.


Cross-correlation peaks in acoustic multi-path conditions

Depending on the circumstances, you can also have other than air borne propagation like through a concrete floor and the support structures of your projector and sensors.  the fastest path isn’t necessarily going to be the direct air path.


Noise removal from audio using FFT/Spectral editing

There exist solutions that are coarsely based on the idea that you propose: Threshold the spectrogram amplitudes to reduce/remove the noise part of the signal. Usually, this will bring about artifacts as e.g. the so-called musical noise and the application determines whether this is acceptable or not. In general, the more the noise is reduced the more negative effects the process has on the quality of the desired signal (i.e. the speech in this case).

As a first attempt I would suggest you to try the Noise Removal algorithm in Audacity (http://web.audacityteam.org/). Probably this will give satisfactory results and is sure worth a try as the software is free. If you're okay with spending money to most probably achieve better results, have a look at RX by iZotope (https://www.izotope.com/en/products/audio-repair/rx/). However, if you know the basics of digital signal processing and want to experiment yourself, a first thresholding algorithm will not be very complex in programming environments like e.g. MATLAB.

If you're interested in the background / signal processing methods that these algorithms are based on I would recommend the books by Vaseghi (http://www.amazon.com/Advanced-Digital-Signal-Processing-Reduction/dp/0470754060/ref=sr_1_1?ie=UTF8&qid=1432211838&sr=8-1&keywords=vaseghi) or Godsill and Rayner (https://www-sigproc.eng.cam.ac.uk/Main/SJGSpringer), the second of which can be downloaded for free from the internet.


Bandpass filter on noisy sinusoid that predicts turning points

Assuming your change of amplitude is very slow when compared to the frequency of the sinusoid, this problem can be posed as finding frequency of a noisy sinusoid. You can take small windows of the signal and then take it's fourier transform. Since we know that the amplitude is almost constant and that there is only one sinusoid, we can look for two peaks in the DFT of the signal with some side lobes(Due to windowing). Your goal would be to find those two peaks, which would reduce the error energy in your estimate. Once you have the frequency and the phase information, finding turning points is very easy.

Hope this helps.


Calculating Beam pattern, Directivity of microphone array

You choose the weights the steer the main beam of the antenna in the particular direction you want - let's call this $\theta_s$ i.e. steering angle. Now, given that you are steering the antenna in a given direction you are also interested in the array's response to incoming signals from other directions. This is the directivity pattern and is a function of the angle $\theta$ of the source signal. For conventional beamforming (non-adaptive) the magnitude of the directivity pattern will reach a maximum when $\theta_s=\theta$.


How to Do Convolution a 2D Signal (Image) and Result of Convolution of Two 1D Filters?

You should show that conv2(f3,x) can be implemented in a separable way. In other words, computing x1 = conv2(f1,x), then applying f2 on the results, x12 = conv2(f2,x1) gives the same result.  Or f12 = conv2(f2,f1) then x12 = conv2(f12,x).

You should be aware of the correct "dimensions" for f1 and f2: transpose of each other, so that one works on "rows" and the other on columns.

This is shown in the following Matlab code:

dataImage = zeros(25,25);
% Create a delta-like image
dataImage(13,13)=1;
% Create a random-kernel image
dataImage(11:13,11:13)=rand(3,3);
imagesc(dataImage);colormap gray

f1 = [0.25 0.5 0.25]'; f2 = [0.25 0.5 0.25];
disp(f1*f2);

f3 = [0.0625 0.125 0.0625; 0.125 0.25 0.125; 0.0625 0.125 0.0625];
disp(f3);

dataImagef1f2 = conv2(f2,conv2(f1,dataImage));

dataImagef1f2 = conv2(conv2(f2,f1),dataImage);

dataImagef3 = conv2(f3,dataImage);

subplot(2,2,1)
imagesc(dataImage);
xlabel('data')

subplot(2,2,2)
imagesc(dataImagef1f2);
xlabel('(data*f1)*f2')

subplot(2,2,3)
imagesc(dataImagef1f2);
xlabel('(f1*f2)*data')

subplot(2,2,4)
imagesc(dataImagef3);
xlabel('f3')

Suggestions for tutorials:

The Matlab Image Processing toolbox guide (900+ pages)
An Introduction to Digital Image Processing with Matlab


Pad image and smooth:

Make image of value 1 same size as original. Pad and smooth:

Divide original image by '1' image:

And crop:

MATLAB demo code:

I = double(imread('lena.tif'));
g = fspecial('gaussian',60,10);
Ipad = padarray(I,[10,10]);
imagesc(Ipad); colormap gray(256); pause;
Ipadg = imfilter(Ipad,g);
imagesc(Ipadg); colormap gray(256); pause;
A = ones(size(I));
Apad = padarray(A,[10,10]);
imagesc(Apad); colormap gray(256); pause;
Apadg = imfilter(Apad,g);
imagesc(Apadg); colormap gray(256); pause;
IA = Ipadg./Apadg;
imagesc(IA); colormap gray(256); pause;
I = IA(11:end-10,11:end-10);
imagesc(I); colormap gray(256); pause;

*Note: * MATLAB imfilter already does zero padding. The pad array step is not needed, it is just for demonstration. Actual code would be:

I = double(imread('lena.tif'));
g = fspecial('gaussian',60,10);
Ig = imfilter(I,g);
A = ones(size(I));
Ag = imfilter(A,g);;
Ismooth = Ig./Ag;


You have a typo in:

result(i) = ( b(2)/a(1)) * audio(i);

That should be:

result(i) = ( b(1)/a(1)) * audio(i);


Behold the power of shortest path.

The solution roughly follows what I did here: Detecting and isolating part of an image.

In this case we first treat the image by limiting every pixel value above the water-level grey to that same intensity.

Then performing shortest path on the inverse edge image created with a difference of Gaussian filter we get:

Resulting in the following segmentation:

My MATLAB functions that are needed for this solution:

http://imageprocessing.com.au/research/code/polarTransform.m http://imageprocessing.com.au/research/code/circularShortestPath.m http://imageprocessing.com.au/research/code/linearShortestPath.m

Code:

clear all
close all

%% Load image
I = imread('CabQx.jpg');
I = double(rgb2gray(I));
figure
imagesc(I);
colormap gray(256);
title('Image');

%% Adjust intensity
Ia = I;
Ia(Ia > 49) = 49;
figure
imagesc(Ia);
colormap gray(256);
title('Adjusted Image');

%% Bandpass Kernel
H1 = fspecial('gaussian',[100,100],8); % *** Changed filter size
H2 = fspecial('gaussian',[100,100],2);  % ***
H = H2-H1;
H = mat2gray(H);
figure;
imagesc(H);
colorbar;
colormap gray(256);
title('Bandpass kernel');

%% Edge Kernel
Hedge = imag(hilbert(H));
figure;
imagesc(Hedge);
colormap gray(256);
title('Edge Kernel');

%% Edge filtered image
G = sqrt(imfilter(Ia,Hedge,'replicate').^2 + imfilter(Ia,Hedge.','replicate').^2);
G = max(G(:)) - G;
imagesc(G);
colormap gray(256);
title('Inverse edge filtered image');

%% Shortest path
% Mid point
[r,c] = size(G);
r = floor(r/2);
c = floor(c/2);
% Min radius and max radius
min_radius = 100;
max_radius = 1000;
% Shortest path
[ path, energy, lastIndex, pathImage ] = ...
    circularShortestPath(G, 100, [c,r], [min_radius,max_radius], [400,720*2]);
figure;
imagesc(I);
hold on;
plot(path(:,1),path(:,2),'r','LineWidth',2);
hold off;
colormap gray(256);
title('Shortest path');

%% Segment
figure
imagesc(I);
hold on;
fill(path(:,1),path(:,2),'r');
hold off;
colormap gray(256);
title('Mask')

%%
J = zeros(size(I));
mask = roipoly(J,path(:,1),path(:,2));
segmentedI = I .* mask;
figure;
imagesc(segmentedI);
colormap gray(256);
title('Segmented image');


Best Approach on How To Sample and Process

As I understand, you want to sample a noisy tone signal at 20dB SNR and in 10Hz-40Hz range with best possible amplitude and frequency accuracy using 12 bit ADC of STM32F4 ARM+DSP chip.

To get the most frequency accuracy you need more samples to be buffered and large FFt sizes for a given samplig rate Fs. This is for frequencies on a given range.

Using high sampling rate (oversampling) may be advantageous for increased SNRQ. See Sigma-Delta converters and Noise shapers for details.

Size of FIR filters: As again it is a direct compromise between computing power and processing quality. STM32F4 is quite powerful chip at 168 Mhz operating frequency and 32 bit FPU.

For FFT, you can write your own code if you have time. I can provide you a butterfly implementation. Or you better use it from some free library.

You should select long fft sizes for better frequency resolution.

please have a look at this document: http://support.ircam.fr/docs/AudioSculpt/3.0/co/Window%20Size.html


I do not see any issues with your approach for the frequency domain computation, given that the following is a circular convolution in the time domain of the two signals a and b:

$$corr(a,b) = ifft(fft(a)fft(b)^*)$$

But to receive each bit you need to evaluate the maximum in each interval to determine if it is positive or negative.

Further in actual application, the received signal would be sampled higher than one sample per bit, with the unknown Doppler offset frequency removed (carrier recovery) and aligned such that your start of the FFT is at the start of the symbol. There are also 20 1023 length symbols per actual data bit in GPS.

In this case the FFT would be rate matched to the received signal. You could also process longer blocks as I show in the example plot below, by simply zero padding the rate-matched prn sequence out to the length of the data block prior to taking the FFT.  In the plot below this was done for 1e6 samples of actual received signals sampled at 20 MSps.


The impulse response characterizes an LTI system completely, in the sense that the response to any input can be computed from the input and from the impulse response  by convolution. That's a consequence of the system being linear and time-invariant.
If we define stability as bounded-input bounded-output (BIBO) stability then stability can be determined from the impulse response. For the system to be stable, the impulse response must be absolutely integrable.
Note that when describing systems by impulse responses, or, equivalently, frequency responses or transfer functions, we focus on the input-output relation. We do not consider the internal dynamical behavior of the system. A state-space model is needed if we not only want to describe the system's input-output behavior but also its internal behavior.


Real time slope and peak detection and calculation

You should start with a bandlimited differentiator (equivalent to a differentiator followed by a low-pass filter). The differentiator will remove the low frequency trend and will respond sharply to your peaks and slopes. The low-pass component will remove noise beyond the cutoff frequency.

You should design your cutoff frequency so that you get clean pulses for your slopes.

Positive slopes will slow as positive pulses; negative slopes as negative pulses, and the peak will correspond to the zero-crossing between positive and negative.

This type of filter is normally implemented as a FIR filter. The number of samples for your filter will then depend on your real-time constraints, the sharpness at the cutoff frequency, and the cutoff frequency itself.


Solution

Impulse train, s(t), of period Ts is used to sample signal x(t):

$$s(t) = \sum_{n=\infty}^{\infty}\delta(t-nT_s)$$

And the sampled signal is:
$$x_s(t) = x(t)s(t)$$

Converting s(t) to frequency domain:
$$F\{s(t)\}= F\{\sum_{n=\infty}^{\infty}\delta(t-nT_s)\}$$

$$S(j\Omega)= \Omega_s \sum_{n=\infty}^{\infty}\delta(\Omega-n\Omega_s)$$

Where:

$$ \Omega_s = \frac{2\pi}{T_s} $$

$$ T_s = \frac{\Omega_s}{2\pi} $$

Now we convert the sampled signal to frequency domain:

$$
\begin{aligned}
F\{x_s(t)\} &= F\{\ x(t)\ s(t)\ \} \\
\\
X_s(j\Omega) &= \frac{1}{2\pi}\ \left[\ X(j\Omega)*S(j\Omega)\ \right] \\
\\
X_s(j\Omega) &= \frac{1}{2\pi}\ \left[X(j\Omega)*\Omega_s \sum_{n=-\infty}^{\infty}\delta(\Omega-n\Omega_s)\right] \\
\\
X_s(j\Omega) &=  \frac{\Omega_s}{2\pi}\ \left[X(j\Omega)* \sum_{n=-\infty}^{\infty}\delta(\Omega-n\Omega_s)\right] \\
\\
X_s(j\Omega) &= \frac{1}{T_s}\ X\left(j\Omega-jn \Omega_s\right) \\
\\
X_s(j\Omega) &= \frac{1}{T_s}\ X\left(j\Omega- \frac{j 2 \pi n}{T_s}   \right) \\
\end{aligned}
$$

Discrete-Time Fourier to Continuous-Time Fourier relationship:

$$x[n] = x(t=T_s n)$$

Take Fourier transform with respect to n to both side

$$ F_n \{ \text{ } x[n] \text{ } \} = F_n\{ \text{ } x(n T_s) \text{ } \} $$
$$X(e^{j\omega}) = \frac{1}{T_s} X\left(\frac{j\Omega}{T_s}\right) $$

Comparing $X_s(j\Omega)$ to $X(e^{j \omega })$, we see that $X_s(j \Omega)$ is a frequency-scaled version of $X(e^{j \omega})$ with the scaling defined by $\Omega=\omega/T_s$. The scaling that makes: $X(e^{j\omega})$ periodic with a period of: $2\pi$ is a consequence of the timing scaling that occurs when $x_s(t)$ is converted to $x[n]$

$$
X(e^{j \omega})=\left.X_s(j\Omega)\right|_{\Omega=\omega/T_s} = \frac{1}{T_s}\ X\left(\frac{j\omega}{T_s}-\frac{jn\omega_s}{T_s}\right)
$$

Properties Needed:

Property 1
$$F\{f_1(t)f_2(t)\} = \frac{1}{2\pi}[ F_1(j\Omega) * F_2(j\Omega)] $$

Property 2

$$F\left\{\sum_{n=-\infty}^{\infty} \delta(t-nT)\right\} = \Omega_s \sum_{n=-\infty}^{\infty}\delta(j\Omega - jn\Omega_s)$$
$$\text{   }Where: \Omega_s=\frac{2\pi}{T}$$

Property 3

$$a\text{ }(f_1(x)*f_2(x)\text{ } )=f_1(x)*af_2(x)=af_1(x)*f_2(x)$$

Property 4

$$f(x)*\delta(x-a) = f(x-a)$$

Property 5: scaling

$$ F\{ f(at) \} = \frac{1}{|a|} F(\frac{j\Omega}{a})$$

Graphs

# From Jupyter Notebook
%matplotlib inline
from numpy import *
from matplotlib.pyplot import *

def FourierSigCT(omega, omega_o=2.5):
    #   this Function "emulates" the spectrum of
    #   a fourier transformed CT signal with max frequency
    #   at +/- omega_o.
    #   Its created by taking a snipplet of a Cosine
    #   between -pi/2 and pi/2, and scaling it to
    #   the cutoff frequency in width..
    cutoff=(pi/2) / omega_o
    if ((omega < omega_o) and (omega > -omega_o)):
        return cos(omega*cutoff)
    else:
        return 0

def FourierSigSampled(omega_o, omega_s):
    t = arange(-10, 10, 0.01)
    n = arange(-5, 5, 1)
    Ts = 2*pi / omega_s
    y = zeros(len(t))
    for j in range(0,len(n)):
        for i in range(0,len(t)):
            y[i] = y[i] + (1/Ts)*FourierSigCT(t[i] - n[j]*omega_s, omega_o=omega_o)
    title("omega_o=%g  omega_s=%g" %(omega_o, omega_s))
    plot(t,y,'.')

t = arange(-3*pi, 3*pi, 0.1)
y = zeros(len(t))

for i in range(0,len(t)):
    y[i] = FourierSigCT(t[i])
title('X(j Omega)')
plot(t,y,'.')

FourierSigSampled(omega_o=2.5, omega_s=2.5*2)
FourierSigSampled(omega_o=2.5, omega_s=2.5*1.5)
FourierSigSampled(omega_o=2.5, omega_s=2.5*3)


Knowing the range of frequencies is a good start and yes a bandwidth filter will certainly help you isolate them but you must consider the overall level of background noise vs signal noise.

Performing an FFT on filtered (bandwidth or otherwise) data will only yield results on the frequencies that are left over! This will either be of use to you or not.

If the number of frequencies that you are after are not that numerous, a collection of IIR (infinite impulse response) filters might be a better choice.

Consider DTMF tones used by telephones and cellphones, typically IIR filters have found their way in software based detection of these dual-tone based signals. There are online resources for calculating each of the  filter's co-efficients. You must use one filter per frequency that you are after.

Here is an iPhone version of source code for DTMF decoding using IIR filters

The file DTMFDecoder.m contains the bulk of the code with the co-efficients defined near the top.

Without knowing any further detail, this is all I am able to help you with.


It seems your needs matches the Dynamic Time Warping algorithm.
You should try it as a metric to compare them.

The idea of changing the size is trickier as you may loose / change data which is important for your comparison.

Maybe what you should do is alter the "Distance" metric in the Dynamic Time Warping algorithm to something which suits your needs.


Can deep learning be applied to video compression?

[EDIT: addition of a March 2017 preprint]

Deep learning already has many applications in video, like enhancement (Deep Convolutional Neural Network for Decompressed Video Enhancement) or semantic analysis.

Recently, there have been some announcements related to video compression, for instance:

Generative Compression (Preprint, March 2017)

  Traditional image and video compression algorithms rely on
  hand-crafted encoder/decoder pairs (codecs) that lack adaptability and
  are agnostic to the data being compressed. Here we describe the
  concept of generative compression, the compression of data using
  generative models, and show its potential to produce more accurate and
  visually pleasing reconstructions at much deeper compression levels
  for both image and video data. We also demonstrate that generative
  compression is orders-of-magnitude more resilient to bit error rates
  (e.g. from noisy wireless channels) than traditional variable-length
  entropy coding schemes.

Video: Deep Learning reconstructing Blade Runner for new video compression techniques
Wyner-Ziv Video Coding using Hadamard Transform and Deep Learning

How important is deep learning here, and what performances are obtained is not clear to me yet.


Reading the amplitude values of a FFT in ImageJ

Looking at the source code, the power spectrum image is using a logarithmic scale, adjusted based on the image's dynamic range such that the minimum value (not less than 1) gives a value of 1 and the maximum value gives a value of 254:

if (min<1.0)
    min = 0f;
else
min = (float)Math.log(min);
max = (float)Math.log(max);
scale = (float)(253.0/(max-min));

...
if (r<1f)
    r = 0f;
else
r = (float)Math.log(r);
ps[base+col] = (byte)(((r-min)*scale+0.5)+1);

Making a broad assumption that at least one of those darker pixel is a null would imply that the computed min is close to 0. With this assumption we can proceed to find max as follows:

$$
\begin{align}
   32 &= \left\lfloor \ln(3.858) \frac{253}{\ln(r_\max)} + 1.5 \right\rfloor \\
   31 &\approx \frac{253 \ln(3.858)}{\ln(r_\max)} \\
   \ln(r_\max) &\approx \frac{253\ln(3.858)}{31} \\
   \ln(r_\max) &\approx 11.02 \\
   r_\max &\approx 61020
\end{align}
$$

Since most of the image is white, I would expect a large spike at the 0 frequency point with an amplitude near $NM = 256^2 = 65536$ (where $N$ and $M$ are the width and height of the image). So $r_\max \approx 61020$ is quite plausible.

So for this image (and this image only, since different images would have a different dynamic range), the approximate scaling would be:

$$
\begin{align}
\mbox{value} &\approx \left\lfloor 22.96 \ln(\sqrt{X_\text{Re}^2 + X_\text{Im}^2}) + 1.5 \right\rfloor \\
  &\approx \left\lfloor 11.48 \ln(X_\text{Re}^2 + X_\text{Im}^2) + 1.5 \right\rfloor
\end{align}
$$


Designing High Pass Filter to Remove DC Offset from a Signal

If the DFT is the Uniform sampling from $ 0 $ to $ 2 \pi $ then the first bin is given by:

$$ x[k] = \sum_{n = 0}^{N - 1} x[n] $$

Namely it is the sum of all the samples.
Hence in order to remove the DC (Mean) all you need is a filter which has zero in its DC bin.
Since, the filtered signal, which is a convolution (Circular) of the DFT of the input signal and the filter will have zero at the first bin which means the sum of the output is zero which means its mean is also zero, as wanted.

Simple and intuitive FIT would be an FIR with the length of signal which removes from the current sample the mean of all samples:

$$ y[m] = x[m] - \frac{1}{N} \sum_{n = 0}^{N - 1} x[n] $$

This is a simple FIR.

Pay attention I assumed you have all the samples.
If you don't, you need to make this FIR casual by only "Calculating" the sum of given samples.


I am trying to construct an aliasing image.

I=rgb2gray(imread('Liver.png')); % a 512*512 image
subplot(1,2,1),imshow(I),title('Original image')  % shown on the left side of figure below

when you implement F=fftshift(fft2(I));, the corresponding axis represent the frequency range of [-255:256]/512; and the cut-off frequency U = V = 256/512 = 0.5;

based on Nyquist Sampling Theorem, the sampling periods in spatial domain should satisfy

delta_x<=1/(2*U)=1;
delta_y<=1/(2*V)=1;

to avoid the aliasing.

So if your bandwidth of recorder is 64 * 64, the corresponding point interval in spatial domain is 1 / (64/256) = 4. Thus,

Ia=imresize(I(1:4:end,1:4:end),[512 512]);
subplot(1,2,2),imshow(Ia),title('Aliased image') % shown on the right side of the figure below

You can observe the aliasing effect at the edge.


How does the RaspberryPi radio hack convert digital GPIO signals to an FM signal?

First, the answer to why you do not see a negative voltage is that the output being digital will range from 0 to the maximum digital voltage at the output (+Vs). This will have a DC offset of +Vs/2 which is simply filtered out with a high pass filter (series cap) resulting in a bipolar waveform with negative voltage after the series cap.

UPDATE: Based on Marcus' link that he provided in his answer, what follows is NOT how they actually did the FM transmitter in the Rasberry PI, but still valid as a possible approach. The actual implementation appears to be much simpler, leveraging an existing PWM spread-spectrum clock controller that is in the Rasberry-Pi. This device is typically used to help pass emission standards by reducing maximum transmitted emissions by spreading them, reducing overall power spectral density. Thus such a device could be used for FM transmission, in that it achieves it's spreading by changing the zero crossings (slightly) of the clock signal being spread. Thus if you can access the control mechanism for how the zero crossings are varied, you can achieve an FM transmitter.

The original response is below to show what could be done if a dedicated spread spectrum PWM clock controller was not available:

The implementation could be done with a numerically controlled oscillator in combination with a sigma-delta DAC. The NCO can provide precise frequency control as well as AM, FM, or PM modulation (you are using FM), and the sigma delta DAC can provide a multi-level analog output from just 1 digital bit.

Since FM is only concerned with zero crossings, it is also feasible that the sigma-delta is completely omitted for a very simple demonstration application in which case the NCO is implemented with a 1 bit output. This would have significant jitter however as the zero crossing will be limited to the crossings of your actual clock vs what is needed in the emulated waveform (as well as high spectral noise requiring more stringent analog filtering). Higher clock rates will reduce this but much more efficient to implement the sigma-delta architecture and in the process minimize that jitter; which is phase noise, which in turn is FM noise; however that jitter is easily quantified/predictable and perhaps with a sufficiently high clock can be shown to not be a concern.

In this post I give more details on the NCO specifically.

Numerically Controlled Oscillator (NCO) for phasor implementation?

As for generating the analog output with modulation, the NCO is combined with a DAC as shown in the figure below which also shows all the modulation knobs (then commonly called a DDS or Direct Digital Synthesizer):

I have not reviewed the Rasberry pi implementation specifically to know if you are making use of an actual n-bit DAC, or if that portion is implemented with a PWM style DAC (Sigma-Delta Modulator), or if the MSB of the phase accumulator is used directly (either of the latter two is what I suspect based on the sketch you made describing the output). In the case of a Sigma Delta DAC, it could be implemented with a pass-band architecture such that the noise shaping is done to maximize the effective number of bits in your FM band (rather than a low-pass approach which would necessitate a very high sampling rate).

Further considerations for DAC reconstruction are given below, not directly related to a Sigma Delta approach but do show how it is possible to generate higher frequency outputs above the sampling rate of the DAC, based on post-filtering after the DAC output.

A very simple Sigma Delta DAC implementation is shown in the figure below. Combining this with the NCO structure above to generate the sample by sample digital word, using the FCW port as the FM modulation control can create an FM signal. As suggested above however, instead of a digital accumlator which is a low pass structure, an implementation can be done with a passband structure that would then allow selective noise shaping to suppress the quantization noise within your FM band of interest:


How to analyze reverb in existing sound recording?

You may be able to reverse engineer the reverb if you have a segment of relatively clean reverb 'tail'. This will allow you to play around and get some metrics such as decay time and frequency response.

Convolution reverb works by recording an impulse response in a particular location/machine after a burst of noise. This is then convolved with the input signal to apply the reverb.


Piano sound generation is complicated by the sound being so familiar to most of us, and the slowly decaying non-stationary non-harmonic nature of overtones meaning that simple attack + looped samples is hard to get right. It is simplified by the low dimensionality of user interaction (88 keys with velocity, 3 pedals) compared to, say, a guitar where there are many dimensions to how the string is picked and how your left hand interacts with the string.
It used to be that digital pianos used synthesis in order to produce tones. My old Roland piano used «SA synthesis»:
http://www.novelmusic.com/roland/Technical/default.htm
As digital memory has become cheaper, my impression is that nowadays, digital piano generation is mostly PCM samples recorded in multiple key (-zones) and multiple strike velocities, with a minimum of resynthesis to smooth things out. I might be wrong.
https://en.m.wikipedia.org/wiki/Inharmonicity


Biquad high-pass filter not working as expected

There is (at least) one error in your routine. It concerns the definition of alpha. The correct formula is

$$\alpha=\frac{\sin(\omega_0)}{2Q}$$

However, you implemented

$$\alpha=\frac{\sin(\omega_0)\cdot Q}{2}$$

But that shouldn't influence the function of the filter as a high pass filter; it just implements the inverse of the given $Q$ value.

Other than that, it's of course important that the variables x1, x2, y1, and y2 retain their value between function calls. I'm not sure if and where you take care of that.

As a final remark, you might want to think about if it's really necessary to recompute the filter coefficients each time a new sample arrives.


Uncertainty Relation between time and frequency

Sampling a signal creates its spectrum replica as Nyquist Shannon theorem stated. As a consequence, to avoid aliasing, one must sample at a rate greater than twice the max frequency $B$ of the signal spectrum. Otherwise, you have aliasing as illustrated in the Figure below.

In your case, $B = 20$kHz hence sampling rate $f_s > 2B$ and the interval between samples $T_s = 1/f_s < 1/(2B)$. You can see that $T_s$ is upperbounded by $1/(2B) = 50$ microseconds.

Note that the longest $T_s$ is not 50 microseconds but $(50-\epsilon)$ microseconds  with $\epsilon$ is arbitrarily small positive quantity.


Difficulties in understanding mutual information concept

What does mutual information (MI) convey?

It indicates that there is a relationship between the two signals- i.e. that they are not independent.  It could be that they are correlated, but the relationship does not have to be that simple.  The higher the mutual information, the stronger the relationship (i.e. the more one of the signals can "tell" you about the other).

  What does this imply and mean in the case of mutual information of error, and in general definition of MI and entropy.

If the error has mutual information with the observation this implies that the receiver did not use all of the information in the observation- i.e. the receiver is not ideal.  If the receiver was ideal it would have used all of the information in the observation, so what was left in the error signal would be independent of the observation.

  What is the interpretation if MI of error is decreasing or getting minimized?

Decreasing MI of error implies that the receiver is doing a better job.  Once the MI is zero it is impossible to do any better.


Can we use ICA (Independent component analysis) on a noisy signal

Yes and no.

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

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

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


Digital audio data that is transmitted wirelessly is usually encoded in some efficient format. It sounds like you're talking about broadcasting raw PCM samples that are encrypted before transmission. That's certainly unusual.

Unless you are already following some kind of standard that I'm not aware of, you're probably best suited by packing your audio data into the payload some higher-level protocol that includes synchronization data to recover framing, the encryption that you desire, and most likely error correction. Your approach of just encrypting raw PCM samples with no error correction has the unpleasant side effect that a single bit error in the encrypted data will turn your entire block of samples into garbage.


Pre-emphasis / De-emphasis filter design

You really can use any filter design scheme that you want. I can't think of any reason why you would need to use an IIR approach, exclusively, so you have a number of options. For FIR filters, you have choices like the frequency sampling and least squares methods (look at fir2 and firls in MATLAB).

For an IIR implementation, you should be able to use the bilinear transform straightforwardly; since the filters you referenced are simple single-pole and single-zero ems, it's probably just as easy to work the math out by hand. One thing to remember is that you'll want to prewarp the frequency axis first to ensure that you get a good approximation across the passband of interest. It might help for you to provide any error that you were getting from MATLAB; the lack of poles shouldn't be problem.

Edit: A few things to address your comments:

I wouldn't say that firls, being a least-squares-based approach, is a step-wise approximation of a smooth curve. Yes, you provide a number of samples of the desired response to the least-squares algorithm, but you should still get a smooth curve out.
I started to work out the math on applying the bilinear transform to the pre-emphasis filter and came across a problem. Its frequenecy response is:
$$
H_{pe}(f) = 1 + j2\pi f \tau_x
$$
Noting that for the frequency response $H_{pe}(f)$, $s = j\omega = j 2 \pi f$:
$$
H_{pe}(s) = 1 + s\tau_x
$$
The bilinear transform maps the $j\omega$ axis in the $s$-plane to the unit circle of the discrete-time approximation's $z$-plane. So, we should be concerned about what happens to the transfer function as $s$ is allowed to vary along the entire imaginary axis. Note that:
$$
\lim_{s \to \infty} H_{pe}(s) = \infty
$$
This means that, for the discrete-time approximation generated by the bilinear transform, as the angular frequency $\omega \to \pm \pi$, the discrete-time filter's magnitude response will blow up to $\infty$, which is probably not what you want.

This problem would present itself for any analog system whose transfer function has a larger-degree numerator than denominator. So it's probably justified why MATLAB issues the following error:

>> [num, den] = bilinear([75e-6 1], 1, 100e3);
Error using bilinear (line 69)
Numerator cannot be higher order than denominator.

(Note that in your example, you passed an empty vector for the denominator; instead, you should pass 1, since the denominator of the transfer function is unity)

So as it turns out, I'm not aware of any analog-mapping techniques that will definitely work for your application. I don't believe the impulse invariant method will work since the pre-emphasis filter isn't bandlimited. You could use one of the FIR design methods, an IIR least-squares technique (e.g. Berchin's frequency-domain least-squares (FDLS) algorithm), or perhaps look into the matched Z-transform method, which I'm not familiar with.


Are two jointly stationary white noise processes independent?

Stationarity (WSS in particular) does not imply independence, even if they are white noises.

That's to say, there exist (at least one) jointly stationary zero mean white random processes which are also dependent. A simple example is this.

Consider a zero mean, wide sense stationary (WSS), white, real random process $x[n]$ with auto-correlation sequence (ACS), $r_{xx}[k] = \sigma_x^2 \delta[k]$, which is passed through a linear time invariant (LTI) filter with real impulse response $h[n] = M \delta[n]$, to produce the output $y[n] = h[n] \star x[n]$ :

$$ y[n] = \sum_k h[k] x[n-k]  = M x[n]$$

Now it can be shown that $y[n]$ will also be zero-mean, WSS white noise with its ACS given by
$$ r_{yy}[k] = h[k] \star h[-k] \star r_{xx}[k] = M^2 \sigma_x^2 \delta[k] $$,

and further, they are jointly WSS and white, as their cross-correlation depends on the lag $k=m-n$;
$$ \mathcal{E}\{x[n] y[m] \}  = r_{yx}[k] = h[k] \star r_{xx}[k] = M \sigma_x^2 \delta[k] $$

However it can be shown that they are not independent, since the necessary condition for their independence does not hold:
$$ \mathcal{E}\{x[n] y[m] \} \neq \mathcal{E}\{x[n] \} \mathcal{E}\{y[m] \} $$

Where the left hand side is the cross-correlation sequence $r_{yx}[k] = M \sigma_x^2 \delta[k] $, as defined above, which is not zero for all $k$, whereas the right hand side is identical to zero for all $n,m$ as both processes are zero-mean WSS.

Hence these two jointly WSS and white random processes $x[n]$ and $y[n]$ cannot be independent and are in fact dependent, in line with intuition as the filtering creates a dependency between them.


Is u[-n] defined from negative infinity to zero as it is folding of u[n] or from negative infinity to -1?

$u[n]$ is the unit-step function :

$$u[n] = \begin{cases} { 1 ~~~,~~~n \geq 0 \\ 0 ~~~,~~~ n < 0 } \end{cases}$$

then

$$u[-n] = \begin{cases} { 1 ~~~,~~~ -n \geq 0 \\ 0 ~~~,~~~ -n < 0 } \end{cases}$$

or stated as
$$u[-n] = \begin{cases} { 1 ~~~,~~~ n \leq 0 \\ 0 ~~~,~~~ n > 0 } \end{cases}$$

Hence, $u[-n]$ is nonzero for the range of $n \in (-\infty,0]$.


How to add a fake thermal effect to images

After experimenting with a couple of libraries, my code finally worked with OpenCV. Code and Output can be seen below.
Code
#this code uses color maps which works as thermal camera filter on the input image

import sys
import numpy as np
import cv2
import matplotlib as mpl
import matplotlib.cm as mtpltcm

def main(argv):

    while (True):
        # read the file
        input_img = cv2.imread('path_to_the_file')

        # convert image into grayscale
        gray = cv2.cvtColor(input_img, cv2.COLOR_BGR2GRAY)

        #initialize the colormap
        colormap = mpl.cm.jet

        #add a normalization
        cNorm = mpl.colors.Normalize(vmin=0, vmax=255)

        #init the mapping
        scalarMap = mtpltcm.ScalarMappable(norm=cNorm, cmap=colormap)

        colors = scalarMap.to_rgba(gray)

        #to display the image
        cv2.imshow("_animal_name",colors)

        if cv2.waitKey(1) & 0xFF == ord('q'):
            break

    cv2.destroyAllWindows()

if __name__ == '__main__':
    sys.exit(main(sys.argv))

Expected Output


Filtering an image with two gabor filters

If you're using just straight 2D convolution, you will get the full output. What I think you want is conv2(A, B, 'same').

To better understand this, 1D convolution between two sequences with length N & M results in a final sequence of length = N + M - 1. So without any truncation your conv2 result will no longer be of size 160x160, but (160+M-1)x(160+M-1).

If you're using Matlab, see 'help conv2', 'help filter2', 'help imfilter' (if you have the image processing toolbox installed). imfilter defaults to 'same'.

I'll add a little more discussion regarding filter sizes when they contain exponential functions.

Given in reality were sampling a function that extends to infinity in all directions, albeit the values become insignificant at some point, the function does extend that far. The implementation you're referencing fixes the limits to 3 standard deviations, because that's good enough for them. I saw another implementation that used 8. If you're scale, or sigma, is very large, you run the possibility of the filter turning out to be larger than the image you'll be filtering against. What do you do then? Because like I stated earlier, straight convolution produces a sequence that is N + M - 1.

You need to grab only the central part of that result. The central part would also be of the same size as the image your filtering. That's what the 'same' flag does.

One other point I'd like to make is that there's nothing stopping you from just creating a filter that is the same size as your image. Instead of using sigma as the size controller, just specify the dimension of the filter. Given convolution is often times done in the frequency domain, you're going to be zero padding and taking a larger FFT anyway. Look for references on fast convolution using the FFT.

Let's say you're image size is 256x256 and you're filter is 32x32. Here's a simple example without using power of 2 FFT's.

N = 256;
M = 32;
K = N + M - 1; % Length of full convolution

x = randn(N, N);
y = ones(M, M);

% Spatial domain convolution (slow)
r = conv2(x, y, 'same');

% Spatial frequency domain convolution (fast)
z = ifft2(fft2(x, K, K).*fft2(y, K, K));
idx = round((K-N)/2+1):round((K+N)/2); % Indices of the central portion
z = z(idx, idx); % Obtaining the central part of the result

subplot(311), imagesc(r), axis('square'), colorbar
subplot(312), imagesc(z), axis('square'), colorbar
subplot(313), imagesc(abs(r-z)), axis('square'), colorbar


Design a filter which passes all frequencies except $\omega=\pm\frac{\pi}{2}$ and plot its pole-zero diagram

The kind of filter you are looking for is a notch filter. Using filter design toolbox in Matlab you can get it as I'm showing you in the following picture:

  This has to be done in z-plane so there must be two poles at +i and −i since they cannot be included in region of convergence. Is my assumption correct?

Yes! This is correct.

Assuming you want to create a causal filter (a right sided sequence with non-zero samples only for n > 0) then your Region of Convergence ROC is :

If you want to understand this better I would recommend you take a look at example 3.1 of Alan V. Oppenheim Discrete-Time Signal Processing 3rd edition.

Normalized frequency response means that instead of dealing with a sampling frequency you make a design independent of it in which you only care about digital frequency. For instance if I say that $w = \pi/2$ then if sampling frequency is 48 kHz (this is $2*\pi$), the maximum frequency that you may get is $\pi$ or 24 kHz and the notch will be at 12 KHz which is $\pi/2$. But on that very same filter if sampling frequency is later 96 kHz then the notch will be at 24 kHz.

The way in which you can do this filter in Matlab is this:

If you select a broader bandwidth (for instance 0.1 instead of 0.001) you will notice how the zero and the pole maintaining the same angle start to separate from each other in the diagram

with sampling frequency 48 kHz:

with sampling frequency 96 kHz:


What do the 1D filters represent when using imfilter?

If I'm not mistaken, a column vector will filter the image across its columns, treating each row independently of the others. Likewise, a row vector will filter across rows, treating all columns the same.

edit: Regarding an example - consider the simple image [1,1,1;0,0,0;-1,-1,-1]. It's constant along its rows (i.e., all the columns are the same) and a gradient along its columns. Let's filter it with a differentiating filter kernel [1,-1] along rows or columns. If we do imfilter(I,[1,-1],'circ') we obtain the zero image: since the image was constant along rows, filtering each row with a differentiating kernel gives the zero image. Oh the other hand, for imfilter(I,[1,-1]','circ') we obtain [1,1,1;1,1,1;-2,-2,-2]: each column gets differentiated independently and since all columns are the same, the resulting image is constant along rows.

As for the 0.25: this is merely a scaling of the whole image. You might as well filter with 1 as a filter weight and divide the result by 4, the effect is the same.


Filtering a signal by modifying the Fourier coefficients and then inverse Fourier-transforming it is a filter design method known as the "frequency-sampling" method.

Many research works on source separation for instance consider such filters because they are convenient (working on spectral characteristics to separate the signals - e.g. the fact that the frequency support of 2 sources is different) easy to implement and still lead to very convincing results (even the infamous "binary masks" :D).

You could compare this method to others, for instance FIR filter design by the window method : compute the ideal IIR impulse response and crop it or window it with a Hann/Hamming/your-preferred-window function. The window method corresponds to multiplying your impulse response by a window, in the time domain. It therefore results in a convolution between the involved Fourier transforms, leading to the well-known Gibbs effect.

For the frequency sampling method, you don't need to find the impulse response, but this comes at the cost that you do not control much of what is happening /between/ the frequency bins: in your example, that would correspond to the frequencies $\frac{\nu}{N} f_s$ where $\nu \notin \{0,1, 2,..., N-1\}$. The actual frequency response of such a filtering process is then an interpolation between the sampled points of the desired frequency response. This issue is not often discussed in the literature, but could lead to big cripples in the resulting frequency response.

There are other methods for designing filters, for instance for designing low/high pass or band-pass filters, but I have less practical experience there. I guess those are anyway not what you are looking for. I mentioned the window method only as a way of comparing what's happening: with the window method, you have trouble with the Gibbs effect, and with the frequency sampling method, you don't really know what's happening between the sampled frequency bins.

This said, for your practical problem, from my own experience, zeroing frequency bins is usually fine. Assuming the way you decide how to enhance this or that frequency bin makes sense (avoiding too many "jumps" in the frequency response, for instance), the result should not exhibit too many artifacts from that process. Note that setting the low values to zeros, instead of keeping the original ones, may however deteriorate what is happening in the frequency region of interest. In source separation as I have experience it, we usually try to minimize some criterion (e.g. mean squared estimation errors), and that can provide us with an optimal way of setting all the values of the Fourier transform. I would argue that this is a better way of dealing with the matter, instead of trying all sort of values for the "missing" ones until the result is acceptable (a method which is most likely to depend on the test signal, hence poorly generalizing to other signals).


How to make sense of audio file's $\textrm{dBFS}$?

However, the dbFS values I'm getting are consistently below -20, i.e., 10% loudness, even at loud parts of the song. Why is this so? Am I calculating the loudness wrongly?

Why should that be a problem? Human perception is roughly logarithmic. Sound that isn't "unbearably loud" shouldn't be much more than 10% of your "acoustic dynamic range".

By the way, with $\frac{}{20}$, this formula gives you amplitude, not power values, so this isn't truly loudness, just amplitude. But since that, in perception, is but a constant factor, these terms can be used interchangably.

  Also, how can I calculate the overall loudness at a point in time? Is it the sum of all the values, or the maximum?

It's the sum, at least approximately: human perception has a different sensitivity for different frequencies, and also, frequencies can mask each other, so the energy in the signal is only proportional to the perceived loudness if we ignore these facts (see: psychoacoustic models).


Why transforming the data to a high-dimensional feature space in which classes are linearly separable leads to overfitting?

@ffriend has a good post about it, but generally speaking, if you transform to a high dimensional feature space and train from there, the learning algorithm is 'forced' to take into account the higher-space features, even though they might have nothing to do with the original data, and offer no predictive qualities.

This means that you are not going to be properly generalizing a learning rule when training.

Take an intuitive example: Suppose you wanted to predict weight from height. You have all this data, corresponding to people's weights and heights. Let us say that very generally, they follow a linear relationship. That is, you can describe weight (W) and height (H) as:

$$
W = mH - b
$$

, where $m$ is the slope of your linear equation, and $b$ is the y-intercept, or in this case, the W-intercept.

Let us say that you are a seasoned biologist, and that you know that the relationship is linear. Your data looks like a scatter plot trending upwards. If you keep the data in the 2-dimensional space, you will fit a line through it. It might not hit all the points, but thats ok - you know that the relationship is linear, and you want a good approximation anyway.

Now lets say that you took this 2-dimensional data and transformed it into higher dimensional space. So instead of only $H$, you also add 5 more dimensions, $H^2$, $H^3$, $H^4$, $H^5$, and $\sqrt{H^2 + H^7}$.

Now you go and find co-efficients of the polynomial to fit this data. That is, you want to find co-efficients $c_i$ for the this polynomial that 'best fits' the data:

$$
W = c_1H + c_2H^2 + c_3H^3 + c_4H^4 + c_5H^5 + c_6\sqrt{H^2+H^7}
$$

If you do that, what kind of line would you get? You would get one that looked a lot like @ffriend 's far right plot. You have overfit the data, because you 'forced' your learning algorithm to take into account higher order polynomials that have nothing to do with anything. Biologically speaking, weight just depends on height linearly. It doesnt depend on $\sqrt{H^2 + H^7}$ or any higher order nonsense.

This is why if you transform the data to higher order dimensions blindly, you run a very bad risk of overfitting, and not generalizing.


Text in image is low or high frequency component?

The edges of the words will probably contribute to a high frequency, the background of the text and the portion within the boundaries of the letters of the text will be otherwise.


How to extract single side of signal spectrum in SIMULINK

If you do not have a single side band separation block-diagram in Simulink, then you would try the following algorithm to separate the upper sideband from the lower sideband given that you have a block called Hilbert transformer.

Given that your original real valued baseband input signal is $x[n]$, then the following signal is called as the analytic signal:

$$  x^+[n] = x[n] + j \hat{x}[n] $$ where $$\hat{x}[n] = \mathcal{H} \{x[n] \} = \frac{1}{\pi n} \star x[n]$$ is called as the Hilbert transform of the signal $x[n]$.

It can be shown that the Fourier spectrum of the analytic signal is:
$$ X^+(\omega) = \begin{cases}
2 X(\omega) &, \text{for } ~~~  0 < \omega < \pi \\
0 &, \text{for } ~~~ -\pi < \omega < 0\\ \end{cases}
$$

Which is 2x the upper side band of the signal $x[n]$


How to crop a region with smooth edge and without holes?

For an image as simple as this, it might work to just perform opening or closing on the image to smooth out the edges before thresholding. Opening will remove any small, disconnected parts of the object, leading to an undersegmented image (smaller than the actual size). Closing will fill in any small holes in the object, leading to an oversegmented image (larger than the actual object size).

However, I think your image here is clean enough that thresholding should still work. Try a smaller value for the thresold, or use a MATLAB function like imbinarize to threshold it for you.

ex. thresholded_image = imbinarize(image, 'global');


Discrete Fourier Transform: What is the DC Term really?

The DC term is the 0 Hz term and is equivalent to the average of all the samples in the window (hence it's always purely real for a real signal). The terminology does indeed come from AC/DC electricity - all the non-zero bins correspond to non-zero frequencies, i.e. "AC components" in an electrical context, whereas the zero bin corresponds to a fixed value, the mean of the signal, or "DC component" in electrical terms.

As far as practical applications go, the DC or 0 Hz term is not particularly useful. In many cases it will be close to zero, as most signal processing applications will tend to filter out any DC component at the analogue level. In cases where you might be interested it can be calculated directly as an average in the usual way, without resorting to a DFT/FFT.


Simple crosstalk cancellation (or signal separation)?

If I understand correctly what you are trying to achieve is usually called Acoustic Echo Cancellation. In many cases, so-called adaptive filtering algorithms, such as the (normalized) least-mean-squares ([N]LMS) algorithm, can be used effectively. In this context, I can really recommend the book Adaptive Signal Processing by Widrow and Stearns. I'm sure you'll find a number of ready-to-use implementations of the (N)LMS algorithm to experiment with on the internet. (Have an eye on the filter order - it needs to be high enough to capture the length of your crosstalk path.)

Furthermore, as you mentioned that you are already working in the frequency domain, you can have a look on algorithms that are called (you guessed it) frequency-domain LMS algorithms. They originate, via so-called block LMS algorithms, from the classical LMS algorithms. In this context, I really liked the paper by Shynk: Frequency-domain and multirate adaptive filtering.


What is the deffinition of convexity defect in image processing?

In simple words, convexity defect is a cavity in an object (blob, contour) segmented out from an image. That means an area that do not belong to the object but located inside of its outer boundary -convex hull (read documentation of cv::convexHull function for details).

The image below shows it better than million of words. Areas between fingers (all marked with arrows) in this schematics of a hand contour are the convexity defects.

You can find the defects using OpenCV's function cv::convexityDefects.


Signal values we will 'miss' between sampling instances during sampling of band limited signals

I don't have a real answer but I have the feeling that this result will help you out: Bernstein's inequality says that, if the signal $x(t)$ is bandlimited to $|f|\leq B$, then $$\left| \frac{\textrm{d}x(t)}{\textrm{d}t}\right|\leq 4\pi B \,\textrm{sup}_{\tau\in\mathbb{R}}|x(\tau)| ,\,\,t\in\mathbb{R}$$ where $\textrm{sup}$ stands for "least upper bound".

I found about this inequality in Amos Lapidoth excellent (and free in PDF format) book "A Foundation in Digital Communication". A proof can be found in M. A. Pinsky, "Introduction to Fourier Analysis and Wavelets".


Protect an IIR filter from being reverse-engineered

Short answer:

You can't. If an attacker can insert a signal that covers the whole bandwidth (e.g. a white signal, or at least one that has no spectral zeros) into the system (and he can do that over an arbitrarily long time, or add up observations), they will get an output, and can through the magic of correlation get the impulse response.


Hidden Markov Models - Identifying Phonemes

You should first appreciate why an HMM would be useful. The problem that you've described is classify a given utterance as either a yes or a no (it would be also useful to have a "neither" class). Because speech is a time-evolving quantity, it needs a classifier that is able to deal with sequences of feature vectors. Different people say "yes" or "no" with different speeds. For example, I might say "no" and the speech signal could be 200ms long or you could say "no" and the speech signal could be 50ms long.

You might want to look into creating windows of a fixed length for your speech signal and extracting certain transform domain coefficients (like FFT, DCT etc.) from it and use these to train an HMM. Alternatively, the mel-frequency cepstrum coefficients (MFCCs) are also very popular in training phoneme based HMMs for speech recognition.


Yes you can think of phases as describing translations.

If you (kind of passively) wrap the phase to ⟨−π,π⟩ you may be introducing discontinuities that are not truly there considering that −π and π are equivalent. Such artificial discontinuities will hinder calculation of things like group delay.

Unwrapping the phase properly is not trivial. One way to do it properly is to factor the system into such small component systems (connected in series) that you can track the phase frequency response of each without any wrapping ambiguities. The phase frequency response of the composite system will be the sum of those of the component systems. Calculated this way, the phase naturally get values outside of ⟨−π,π⟩. Computationally, the factorization can be done by finding the z-plane zeros of the system. My guess is that numpy.unwrap does nothing of the sort but simply looks for big jumps of the wrapped phase as function of frequency and tries to make those jumps shorter by adding or subtracting 2π. That will easily fail do do a proper job if the frequency response has two zeros at the same frequency, as that gives a jump of 2π that phase unwrapping ought to preserve.


How to calculate signal which is not changed by a filter?

If a filter with transfer function $H(f)$ passes a signal $x(t)$ with Fourier transform $X(f)$ unchanged, then $H(f)X(f) = X(f)$ for all $f$, and
so $H(f)$ must have the property that
$$H(f) = 1 ~ \text{for all frequencies}~ f~ \text{for which}~ X(f) \neq 0.$$
Thus, an arbitrary filter will not enjoy this property for very many
signals, perhaps for no signal at all.  Indeed, if $H(f_0) = H(-f_0) = 1$ for
some $f_0$, then the filter will pass the sinusoid
$A\cos(2\pi f_0t + \theta), -\infty < t < \infty$ unchanged, but it will
not pass a sinusoidal pulse such as
$$x(t) = \begin{cases} A\cos(2\pi f_0t + \theta),
& -T < t < T,\\ 0, &\text{otherwise,}\end{cases}$$
without distorting it. Of course, it is possible that $H(f) \neq 1$ for
all $f$ and so even passing a monochromatic sinusoid is not possible
(cf. @hotpaw's comment).
In short, the subspace referred to by @AtulIngle might well be
the trivial subspace of zero dimension containing only the
zero signal.
It is also worth noting that if $H(f)$ has constant nonzero value
over an interval, then the filter is physically unrealizable for
essentially the same reason that an ideal low-pass filter is
unrealizable.

So, the OP's question is likely to have no answer at all because the
signals $S^\prime$ that he seeks to find may not exist at all for
the given filter. Perhaps something else was intended to be asked?


Bilinear Interpolation on A Subsection Of An Image

Bilinear interpolation can be viewed as traditional sample rate conversion (insert 0s between existing samples, convolve with a triangular kernel in each dimension). Somewhat unusually, the size of the (continous-time) triangular filter kernel is not scaled relative to resampling ratio before it is sampled; each sample is always the result of 4 neighbouring input sample (unless it falls exactly on one of them). This makes it quirky seen with dsp goggles but maps very efficiently to compute hw.

When doing convolution you need overlap between separate chunks, or you will get some kind of discontinuity. For bilinear interpolation the overlap need only be one sample wide.


Please check the following for Iris Normalization:

  Nithyanandam. S, Gayathri K.S., Priyadarshini P.L.K. - A New IRIS Normalization Process For The Recognition System With Cryptographic Techniques

  Arvacheh E. - A Study of Segmentation and Normalization for Iris
  Recognition
  Systems


How to Use Maximum a Posteriori Probability (MAP) in Classification Task

I will try to give you some intuition into it by a different example.

Think we have 3 machines which can generate the numbers 1, 2, 3.
The first machine generates the number 1 with 80% and the numbers 2, 3 with 10% each.
The second machine generates the number 2 with 80% and the numbers 1, 3 with 10% each.
The third machine generates the number 3 with 80% and the numbers 1, 2 with 10% each.

One of the machines is chosen, you don't know which and the generated number is 2.
What machine would you bet it would?

The question above is the likelihood with no prior knowledge.
Hence, given the number 2 the most likely is the second machine.

Yet,
What happens if you are told that the chosen machine would be chosen following this rules - 98% the first machine, 1% the second and 1% the third.

Now, what are the chances?

Seeing 2 from machine 1 are 0.1 * 0.98 = 0.098.
Seeing 2 from machine 2 are 0.8 * 0.01 = 0.008.
Seeing 2 from machine 3 are 0.1 * 0.01 = 0.001.

Now a better choice would be machine 1.
This is the MAP, because we took into account both the likelihood and the prior knowledge.


As with any question like this, the answer is: It depends.

What does it depend on? Your signal model.

If your signal model generates $X$ and $Y$ axis velocities independently from each other so that there is no transfer function between the two channels, then two 1D Kalman filters will work the same as one 2D Kalman filter.

If your signal model allows interaction between $X$ and $Y$ axis values, then you lose that captured interaction between the axes, which can degrade the performance of the two 1D Kalman filters over a single more accurate 2D Kalman filter.

To give any more detailed answer will require you to state your signal model.

The edit really doesn't give me enough information, but let me try to update my answer.

I'm assuming, from what you say, that the latitude and longitude are estimated from a nadir image.  If that's true, then unless the image pixels are precisely aligned with the lines of latitude and longitude, you will get cross-coupling between your latitude and longitude estimates.

That cross-coupling means you'll be better off using a 2D Kalman filter.


Is document image binarization a closed research field

First, in science, a field is rarely closed, sometimes asleep only. Resistance to low-contrast, real-time, badly  scanned, composite documents/writers or from aging medium seem to remain challenges, in similarity with other digital data: robustness, speed, poor acquisition, source separation or error concealment (the same topics, with generic names) are common for other data source.

As many data fields, this domain is being pervaded by (deep) learning techniques threatening the old handcraft (eg Learning Document Image Binarization from Data, Yue Wu et al.,  2015).

DIBCO 2017, the ICDAR 2017 Document Image Binarization COmpetition, is still on, a good sign altogether, and a source of challenges.7

As a personal note, I have been very interested in the analysis of old written documents, and their authentification, based on:

binarization from alternative light sources (multispectral, X-ray, ultraviolet), especially to detect previously scrapped, washed off  or erased writings, like from palimpsests (Binarization of MultiSpectral Document Images, 2015)
accurate analysis of non-written elements like the texture of paper fabric


basic hsb skin detection, neon illumination

This actually might work better using a simple generative model in RGB rather than HSV.

Get a training image or several training images with some skin.
Manually select the skin pixels (e. g. by creating a binary mask)
Compute the mean and covariance of the skin tone in RGB (each should be 3-element vectors)
For an unknown pixel, compute its Mahalanobis distance from the mean, using the covariance.
Classify it as skin if the distance is less than a threshold.
Tune the threshold for best performance.

Edit:
I don't know if OpenCV has a function to compute covariance, but I can tell you how to do it yourself. Let's say you have $n$ RGB pixels.  You put them into an $n$ x 3 matrix, let's call it $P$.  Then compute $m$, which is the mean RGB vector by calculating the average of the columns.  $m$ will be a 1 x 3 vector. Subtract $m$ from every row of $P$ and call the resulting matrix $Q$.  Now to compute covariance, all you have to do is multiply $Q$ by the transpose of itself: $C = Q'Q$. Make sure that $C$ is 3 x 3.

Edit2:
The values you are getting seem to be too large.  To get the maximum covariance create the following matrix:

255 255 255
 0   0   0

and compute covariance of that.  You should get a matrix where every value is approximately
32513.  So make sure that your pixel values range from 0 to 255, and make sure you copy them into floats or doubles correctly.  The Mahalanobis distance is in the units of variance, so for that the numbers should be small.  Your threshold for skin classification should probably be less than 4.


Missing f0 is a common situation. Pitch could be defined the lowest common multiple of the spacing between spectral peaks, and they seem to be spaced by 82 Hz or 164 Hz here - so there's nothing wrong with your signal and spectrum. Harmonic sum or product should deal with this case. Could you please post plots of your harmonic sum or product? Which result do they give you?

You can also try looking at another pitch detection method like the Average Magnitude Difference Function or YIN.


According to their model given by Eq. (2), voiced speech is assumed to consist of a weighted sum of sinusoids with frequencies that are multiples of the fundamental frequency. They assume that in each frequency band $k$ there is one dominant harmonic with normalized angular frequency $\Omega_h^k$, which is the harmonic that is closest to the frequency corresponding to the $k^{th}$ frequency bin of the FFT. The frequency of this harmonic is used to estimate the current phase in band $k$ from the previous phase value in that band according to Eq. (3).


Why is the concept of a "state covariance matrix" necessary in estimation?

A resistance isn’t a particularly dynamic state, should be an unknown constant, but you might have a bin of resistors where they might vary. Taking one “randomly” out the bin makes it a random variable.  The bin of resistors will have a mean value, so perhaps that mean constitutes a state variable.  Perhaps someone starts putting resistors in the bin from a different source where they tend to run high, so the mean changes.

You don’t “have” to assume that the state is random if the model isn’t really random. Later you will find out in your class that states do not have to be stochastic. The KF equations will work for deterministic states. The measurement does need to have noise. If it doesn’t have measurement noise, you use a state observer, not a KF.

A random state is actually very general and modeling errors are often successfully considered noise.

Typically we use states because they encapsulate all the previous states. You use the current state to move to the next state. you don’t need to consider any past states to go to the next state. Some people will say a state has memory. If the state is stochastic, it has the Markov property. Your measurement depends only on the current state.

The idea of randomness has many interpretations. It is part of what makes the topic interesting.