Signal Processing/FFT gives very high magnitudes for low frequencies

A very simple solution to your problem is to use a window on your dataset prior to taking the FFT. This classic paper by fred harris details the scalloping loss for all the popular windows. Note the following from this paper:

The scalloping loss for a rectangular window (meaning no window, which I believe is what you are doing) is 3.92 dB.   The minimum listed in Table 1 in his paper is 0.83 dB using a "Minimum 4-sample Blackman Harris". This presents a very simple solution in that you multiply your time domain by the appropriate window function prior to taking the fft, reducing the amplitude variation from 3.92 dB (36.4% drop from peak at bin centers) down to 0.83 dB (9.12% drop)!

Further you can combine this with overlap-add techniques if you are block processing streaming data, which is also further explained in the paper.

https://www.utdallas.edu/~cpb021000/EE%204361/Great%20DSP%20Papers/Harris%20on%20Windows.pdf

Also if your requirements is simply for graphical reasons for a swept sinusoidal signal, do know that the signal is not actually "lost" but distributed to adjacent bins. You can recover most of the signal by root sum-squaring the adjacent bins (or all the signal by rss'ing all the bins and displaying this at the max signal location).  Doing this, for a sinusoidal signal, will eliminate all inter-bin variation.

Alternatively similar to the previous paragraph for a streaming application where we want to compensate the output for "presentation" purposes on a graph, predistortion (equalization) of the passband ripple prior to the fft would be simple to do, as long as you have more data available than the size of your FFT. A 3 tap FIR compensator would bring the ripple from 3.9 dB to less than 1 dB. This technique could also be combined with the windowing mentioned to bring the ripple down well below the 0.83 dB provided by the window.  Without windowing, each bin of the FFT output is a Sinc filter at the bin's center frequency with the nulls at all the other bins. Recognizing this, we could use an inverse Sinc approach as a compensator, following the design for a 3 tap inverse Sinc compensator given as my response to this post:

how to make CIC compensation filter

Trying out a first pass of this approach on a simple 4 pt DFT (to match the plots above), I was able to achieve the result shown in the plot below showing the scalloping loss of -3.92 dB as listed in fred harris' paper, and then the reduced deviation after using a simple 3 tap FIR prefilter.  This can be reduced further with an IIR approach of similar order, or more taps but at that point I would likely use a windowing approach (as the span of the taps is beyond the FFT data length).


Given a impulse responce can we derive a state space model?

you might need some procedure like the prony method to convert your impulse response (as a bandlimited and uniformly-sampled signal) into a sum of decaying exponentials and from that obtain a sum of partial fractions of first and second-order transfer functions, and from that derive poles and zeros.  from that can be derived a state-space model.


Using Doppler velocimetry with digital signals such as Bluetooth

If your question is whether you could ever use the Bluetooth radio in an iPhone to perform doppler velocimetry, I'm going to go out on a limb here and say no.

Most mobile devices are incredibly  modular. This means that the BT/WiFi radio (they are often combined into a single chip these days) is designed to handle all of the RF signal processing necessary, and then hand over only the digital data received, but never the analog signals.  By the time the signal is turned into digital data, any effects of multipath, channel noise, doppler shift, and other analog phenomena has been removed.  It is highly unlikely that any mass-produced cell phone would have access to the analog signal being received by the individual radios.


Minimum number of Poles and zero of transfer function H(z)?

Restricting the discussion to LTI systems with rational transfer functions, there are various conditions on the placement of poles and zeros  of $H(z)$ implied by several properties of the system under concern. Using some of them helps in this question.

Prop 1: Total number of poles equals total number of zeros (including those at $z=0$ and $|z|=\infty$) ,i.e., There are equal number of poles and zeros for any $H(z)$.

Rule 2: For real systems (which has real impulse response $h[n]$) poles (or zeros) exist as either purely real or as complex conjugate pairs.

Now applying these two properties to your system which is real and causal (how?) we reach the following conclusion:

1- If ROC of $(1-\frac{1}{2} z^{-1}) H(z)$ is the entire $\mathcal{Z}$-plane (exc $z=0$), then $H(z)$ has a pole at $z=0.5$

2- If $x[n]=(1+j)^n$ produces zero output, then $z=1 \pm j$ are two zeros of $H(z)$

Are there any more poles of $H(z)$ other than those at $z=0$ ? No. Because if there were so, then part one would indicate it in the ROC of $(1-0.5 z^{-1})H(z)$ which is given as entire $\mathcal{Z}$ plane indicating that there are no more poles of finite, nonzeros values.

Are there any more zeros of $H(z)$ ? The given information is not sufficient to conclude that. There could be more zeros. If they were real then they could be even or odd numbered, if they were complex then they must be in complex conjugate pairs. We only know one complex pair of zeros at $z=1 \pm j$. Note that for each added zero to $H(z)$ there would be one more added pole at $z=0$ (in addition to the existing pole at $z=0$)

Hence we conclude: $H(z)$ has two poles at $z=0$ and $z=0.5$ and at least two zeros at $z=1+j$ and $z=1-j$. If $H(z)$ had any more zeros at $z_k=\alpha_k$ then there would be equal number of poles at $z=0$. And if the zeros were complex they would be in complex-conjugate pairs (and the added poles would be doubled)


Estimate "speed" of repeating signal pattern

You may be able to use a modified transform to compute change in a time stretch parameter (I will use $\zeta$) versus time (t). Your success in being able to do this depends on the cross correlation properties of your pattern with time stretched versions of itself. Let me explain:

First notice that the Fourier Transform is a correlation over all possible frequencies (here using the simpler form with complex tones $e^{j\omega t}$, once you know that $e^{j\omega t}$ represents a single tone in frequency this form is so much cleaner than expressing it with sines and cosines; since a cosine is two complex tones: $2cos(\omega t) = e^{j\omega t} + e^{-j\omega t}$).

$$F(j\omega) = \int_{t=-\infty}^\infty x(t) e^{-j\omega t}dt$$

Observe that in this process, you set $\omega$ to every possible value, and for each $\omega$ you compute the correlation of your waveform with that frequency, where what I refer to as "correlation" generally is a complex conjugate multiplication and integration process (and digitally with the DFT we do a multiply and accumulate over a finite length):

$$Corr = \int x(t) y^*(t) $$

Similarly, the short term FT (STFT) we introduce a time variable $\tau$ where we sweep a window function (which selects a portion of our waveform) to access the change in frequency versus $\tau$:

$$F(\tau, j\omega) = \int_{t=-\infty}^\infty x(t)w(t-\tau) e^{-j\omega t}dt$$

Where $w(t)$ is a window function, and we have the trade space of time resolution and frequency resolution to contend with. (The simplest window is a "rectangular window" which is all one over a specified duration of time).  If our signal is changing very slowly with time, we can use longer windows in the STFT to achieve fine frequency resolution. but if it is changing quickly, and we want to observe that then we must use shorter windows resulting in courser frequency resolution. Within the choice of window itself we have the trade space of resolution versus sensitivity (dynamic range): a rectangular window has the finest frequency resolution but very poor sensitivity: the sidebands roll-off at 1/f. Higher performance windows (Kaiser, Hamming, Blackman etc) can have much better sideband rejection (better sensitivity) but always at the expense of the main lobe width (resolution).

To the degree that your base pattern, what I will call m(t), is sufficiently uncorrelated to time stretched versions of itself ($m(t/\zeta)$), you can modify this transform to be what you are looking for. "Sufficiently" here defines the sensitivity range of your approach. Assuming it is sensitive enough for the actual time stretching you are experiencing, you can then use this appraoch to determine how your actual signal $x(t)$ correlates to to different time stretching of your base pattern. You will have the same trade space as in the STFT: in this case $\zeta$ resolution versus time resolution. If the $\zeta$ is changing very slowing then you can integrate over a long duration and get a fine resolution for $\zeta$ up to the maximum cross correlation properties your base function allows, but if changing quickly and you want to measure that fast change then you will need to use shorter integration lengths.

That said if you want to go down this path I suggest:

Step 1) Analyze the cross correlation property of your signal $m(t)$ with time stretched versions of itself using:

$$rr_{t \zeta} = \int_0^T m^*(t)m(t/\zeta)d(t)$$

Where $\zeta$ is the parameter stretching your time, when $\zeta=1$ you get the autocorrelation at time=0, for values of $\zeta$ less than 1, your waveform $m(t/\zeta)$ will be stretched in time relative to $m(t)$, and for values of $\zeta$ greater than 1, the waveform will be compressed in time. If your waveform is real then you can ignore the complex conjugate, however this is critical to include for complex waveforms, so I included it in the formula.

The amount of cross correlation will indicate the usability of this approach and your measurement sensitivity and precision. Check this with increasing values for T (integration time) to access usability range for T (which would then provide guidance on window lengths). Assuming the cross-correlation is low enough for your use, proceed to step 2.

Step 2) Use the modified transform to determine $\zeta$ versus time for your actual capture waveform:

$$F(\tau, \zeta) = \int_{t=-\infty}^{\infty} x(t)w(t-\tau)m^*(t/\zeta)dt$$

Where:
$x(t)$ is your received signal

$w(t)$ is a window function (rectangular, Kaiser, Hamming etc)

$\tau$ is your time parameter for the rate of change in the stretching of your base pattern

$\zeta$ is the time stretch parameter for your base pattern.

The result is a 3 dimensional plot with $\zeta$ on the horizontal axis, $\tau$ on the vertical axis, and the magnitude result for your correlation along the z axis. Where the magnitude is maximum indicates the best estimate for $\zeta$ at that time $\tau$.


Can the Walsh Hadamard transform be computed for non power of 2 sizes?

The Walsh-Hadamard transform requires a Hadamard matrix -- which has entries $\pm 1$ and
whose rows are orthogonal vectors.  A $n \times n$ Hadamard matrix is known to exist
for $n = 2$. Larger Hadamard matrices can exist only if $n$ is a multiple of $4$, though
it is not known if there is a Hadamard matrix for every multiple of $4$. However,
there is a recursive construction that can be used to construct a $2^m \times 2^m$
Hadamard matrix from a $2^{m-1}\times 2^{m-1}$ Hadamard matrix, and this structure
allows for the use of the Fast Hadamard Transform algorithm, which reduces the
computational cost from $2^{2m}$ additions and subtractions to $m2^m$ additions and
subtractions (just like the $N$-FFT reduces the number of multiplications from
$O(N^2)$ to $O(N\log N)$.

With this as background, the answer is Yes, you can use a Walsh-Hadamard transform
of length $N$ if you can find a $N\times N$ Hadamard matrix, but your choices
for $N$ are necessarily restricted. Also, fast algorithms may not exist
for your choice of $N$, though some speed-up is usually possible.

Note that Walsh-Hadamard transforms (WHTs) do not support cyclic convolutions, but
they do support what is sometimes called Poisson convolution. If $H$ denotes
a $2^m\times 2^m$
Hadamard matrix and $\mathbf x$ and $\mathbf y$ are vectors of length $2^m$ with \WHTs $\mathbf xH$ and $\mathbf yH$,
then the inverse WHT of the term-by-term multiplication of the entries
in $\mathbf xH$ and $\mathbf yH$ can be described as follows.  Suppose that the entries in
$\mathbf x$ and $\mathbf y$ etc are indexed not by integers $0$ through $2^m-1$ but rather
by the $m$-bit representations of these numbers.  Thus, we talk not of
$x[k]$ as the $k$-th entry in $\mathbf x$ but rather of $x[\mathbf k]$ where
$\mathbf k$ is the $m$-bit representation of $k$.  Then, the iWHT of the term-by-term
product of the entries in $\mathbf xH$ and $\mathbf yH$ has $\mathbf k$-th entry
$$\sum_{\mathbf i} x[\mathbf i]y[\mathbf k\oplus \mathbf i]$$
which is eerily reminiscent of
$$\sum_i x[i]y[k-i]$$
and even more so if one notes that modulo two, addition and subtraction
are the same and so that $\mathbf k\oplus \mathbf i$ could as well
have been written as $\mathbf k\ominus \mathbf i$.


Why level of noise can be magnified twice through each numerical differentiation?

Simply put: take an independent identically distributed  Gaussian noise (one observation in blue, left) and its gradient  (in green, right). The amplitude is roughly multiplied by the norm of the gradient filter (here $h=[1\,-1]$), which is $\sqrt{2}$, which you can see from the picture.

Their average energy  $E$ thus differ by a factor of $|h\|^2=\sqrt{2}^2 = 2$. As said by Olli, a similar question was answered in 2nd order edge detectors more susceptible to noise?. The part "can be magnified twice" is important, since this does not happen for all noises or all difference filters.

More generally, a random signal $x(t)$  passing  through a time-invariant filter with frequency response $H(f)$ produces a signal $y(t)$ with  PSD:

$$ S_y(f) = S_x(f) |H(f)|^2\,. $$

Since the Fourier domain is energy-preserving, you get that

$$\|H\|^2\ =|h\|^2=2\,.$$

This is a reason why people have developed noise-robust differentiation, often combined with smoothing (eg Savitzky-Golay filters). See examples at
Differentiation.

Note that in the continuous case, you should take care of non integrable processes, see for instance Variance of white Gaussian noise.


What that statement is saying that

a. Correlation is performed the same way as one would perform convolution (you must implicitly know that one of the sequences is conjugated and time reversed to express a correlation as a convolution, as it was not stated there),

and

b. Mathematically the convolution is performed using this relationship

$$x[n]*h[n] = \mathscr{F}^{-1}\left\{\mathscr{F}\left\{x[n]\right\}\cdot\mathscr{F}\left\{h[n]\right\}\right\}$$

Where $*$ denotes convolution, $\mathscr{F}\left\{f[n]\right\}$ denotes the Discrete Fourier Transform of $f[n]$, and $\mathscr{F}^{-1}\left\{F[k]\right\}$ denotes the Inverse Discrete Fourier Transform of $F[k]$.

Note that these operations with Discrete Fourier Transforms actually yield the "circular" convolution of $x[n]$ and $h[n]$, which is equal to the linear convolution of $x[n]$ and $h[n]$ only under a certain condition: there must be enough zero padding to avoid aliasing.


Why does DTFT start and end at the same magnitude?

DTFT is a periodic waveform;

$$H(e^{j \omega}) = H(e^{j (\omega + 2 \pi k)}).$$

Hence for every frequency interval $[\omega_1, \omega_2]$ of length $2\pi$ it will repeat itself; i.e., beginning and ending at the same magnitude (and phase) on that interval.


how to set Equalizer's coefficient using generalized eigenvector.

The generalized eigenvalue problem is given by

$$Bw=\lambda Cw\tag{1}$$

where $\lambda$ is the generalized eigenvalue of the matrices $B$ and $C$. Multiplying $(1)$ from the left with $w^H$ (with $^H$ denoting the Hermitian conjugate) and dividing both sides by $w^HCw$ (assuming that this term is non-zero), we obtain

$$\frac{w^HBw}{w^HCw}=J(w)=\lambda\tag{2}$$

This shows that the value of the objective function $J(w)$ (which is a generalized Rayleigh quotient) equals the generalized eigenvalue $\lambda$. Consequently, $J(w)$ is maximized by the eigenvector corresponding to the maximum eigenvalue, and it is minimized by the eigenvector corresponding to the minimum eigenvalue.


How to compare multiple polynomials in MATLAB?

Based on what you share, this is the best answer I can give.
I see two methods:

Method #1,

Instead of making N polynomials (if N is the number of ribs), you could fit a curve for all the ribs at once. Of course, for that, you will need to have a common (0,0) coordinate, in other word, stack them one over the other on your image. Since you are able to make a distinction between the ribs, that means that you can already isolate them, so that shouldn't be an issue.

Or if you prefer, you can generate a new cloud of points with the polynomials you already found.

coeffs = [0.0010126, -1.5981, 1074.8; 0.0024871, -3.3327, 1866.3];
x = linspace(0,1000,1000); % fictives points
merged_x = [x x];
merged_y = [ polyval(coeffs(1,:), x), polyval(coeffs(2,:),x)];
coeffs_final = polyfit(merged_x, merged_y, 2);
figure
hold on
plot(polyval(coeffs(1,:), x), 'red')
plot(polyval(coeffs(2,:), x), 'blue')
plot(polyval(coeffs_final, x), 'black')

Method #2

You could fit 3 polynomials to find the values of the coefficient you already have found. If we take the first polynomial you provided, then
$a_2=0.0010126$ $a_1=-1.5981$ $a_0=1074.8$

You could fit a curve that would follow all $a_2$, another for all $a_1$, and same for all $a_0$. This would yield 9 new coefficients (3 polynomials) that will allow you to find your polynomials based on the rib position.

Here's an example

coeffs = [ 0.0010126     -1.5981       1074.8;
    0.0024871     -3.3327       1866.3;
 0.002662      -3.7067       1631.4;
 0.0025574     -4.5129       2148.4;
 0.002114      -3.3457       2195.6;
 0.0015383     -2.6023       1682.9 ];

x = linspace(0,1000,1000); %fictive value
ribs_count = size(coeffs, 1);
a2 = polyfit(1:ribs_count, coeffs(:,1)', 2); % careful to the ' symbol that transpose the column into a row
a1 = polyfit(1:ribs_count, coeffs(:,2)', 2);
a0 = polyfit(1:ribs_count, coeffs(:,3)', 2);

figure
for rib=1:ribs_count
    subplot(6,1,rib);
    hold on
    new_coeffs = [polyval(a2,rib), polyval(a1,rib), polyval(a0,rib)];
    plot(x, polyval(new_coeffs, x), 'red')
    plot(x, polyval(coeffs(rib,:), x), 'blue')
end


Why calculate negative frequencies of DFT?

Mainly because its easier. The FFT is a specific algorithm to calculate the DFT. However, it only works if you calculate ALL frequencies (regardless if you want them or not). It takes in N complex values and it spits out N complex values. So the FFT can be used to evaluate the your first equation but not your second.

In your example, that's a trivial difference; but for larger FFT sizes this makes a substantial difference. It's a lot faster to calculate the values at negative frequencies with FFT and simply throw them away then calculating only the values at positive frequencies with a different algorithm.

This being said, there are ways for making use of the redundancy:

You can pack two real signals into a single complex, do a single complex FFT and split this into odd and even parts to get the DFTs of both input signals
There is a way to do a N point real DFT using N/2 point complex FFT. It just needs one extra stage of "unwinding" the result.


Can a (vanilla) Kalman Filter's Observation Matrix $H_k$ depend on the state vector $x_k$?

The non-linear system model that the unscented and extended Kalman filters work with (omitting control) is the following:
$$x_{k+1} = F(k, x_k, v_k) \\ y_k = H(k, x_k, n_k)$$
where $F$ is the state transition function and $H$ is the observation function.  The linear Kalman filter is just the special case where $F$ and $H$ have the form:
$$ F(k, x, v) = A_k x + v \\ H(k, x, n) = B_k x + n$$ for sequences of matrices $A_k$ and $B_k$.  In particular, the functions $F$ and $H$ are linear.

The following is a bit of a pedantic distinction (particularly in the non-linear case) but while the expression $F(k, x_k, v_k)$ clearly depends on $x_k$, the function $F$ itself does not.  In the non-linear case this is no limitation.  In the linear case, if the matrix $B_k$ depended $x_k$ then the function $H$ would no longer be linear.


DC component in accelerometer data - filter before or after XYZ Euclidean sum

I notice that after the butterworth filtering the high frequency components are attenuated. That suggests that somehow you are actually low pass filtering. There might be a mistake in your code.
A sliding window averaging is a crude low pass filter. Therefore, given that you have both a low pass filter (window averaging) and a high pass filter, you might try combining both operations by using a band pass filter.
You have to remove the DC offset from each channel before taking the magnitude (euclidean sum).  Consider a simple example - the offsets are $[1,2,3]$ and the actual acceleration vector is $[1,1,0]$. The combined magnitude is $\sqrt{(1+3)^2 + (1+2)^2 + (0+3)^2} = \sqrt{34}$. But if the values for two axes were flipped, say $[1,0,1]$, you would get  $\sqrt{(1+3)^2 + (0+2)^2 + (3+1)^2} = \sqrt{36}$. Ideally the magnitude should be invariant to rotation of the axes.


What is the effect of aliasing on the magnitude of the autocorrelation?

Decimating before calculating the autocorrelation, in the presence of noise, is inferior to calculating the autocorrelation using the full dataset. Assume that the signal of interest is embedded in white noise. The vector $x[n], n = 0, 1, ..., N-1$ consists of samples from a discrete random process. The autocorrelation function of the vector $x[n]$ is:

$$
A_x[k] = \frac{1}{N-k} \sum_{i=0}^{N-1-k} x[i]x[i+k]
$$

That is, $k$ is the lag used for the autocorrelation calculation. In your proposed scenario, you are decimating the autocorrelation function output by a factor $D$ (i.e. you are only calculating the function for lags $0, D, 2D, ...$) and comparing that result to the autocorrelation function of $x[n]$ decimated by the same factor $D$. Let $x_d[n]$ be the decimated sequence; its autocorrelation function is:

$$
A_{x_d}[k] = \frac{D}{N-k} \sum_{i=0}^{\frac{N-1-k}{D}} x[iD]x[(i+k)D]
$$

(for simplicity here, I have assumed that $D$ is a factor of $N$ in the above equation)

Your inquiry can be written as:

$$
A_x[kD] \stackrel{?}{\approx} A_{x_d}[k]
$$

$$
\frac{1}{N-kD} \sum_{i=0}^{N-1-kD} x[i]x[i+kD] \stackrel{?}{\approx} \frac{D}{N-k} \sum_{i=0}^{\frac{N-1-k}{D}} x[iD]x[(i+k)D]
$$

Looking at this qualitatively, the summation on the left hand side has more terms than its counterpart on the right side. If $x[n]$ is second-order stationary, then the expected value of each term in each sum is the same; the act of averaging multiple samples that have the same expected value increases the signal to noise ratio. Stated a little differently, you can think the terms in each sum as samples from a new random process:

$$
y[n] = x[n]x[n+kD]
$$

Since the noise present in $x[n]$ is white, the expected value of $y[n]$ is the true autocorrelation of the signal of interest at lag $kD$. Therefore, we would like to accurately estimate the expected value of $y[n]$. Our method for doing so is by calculating a sample average; it can easily be shown that the variance in the sample average estimator decreases given a larger sample size, converging to the actual expected value as the number of samples tends to $\infty$.

So, if there is white noise present in the signal (which is often the case), you're going to get a better estimate of the underlying signal's second-order statistics by using a larger sample size in the calculation (this might sound intuitively obvious). In the context of your two approaches, this is accomplished by using the full, non-decimated signal in the autocorrelation calculation and decimating afterward (i.e. only calculating the result for certain lag values).


Are there equations for placing sounds in the stereo field using delays?

Spherical head, distant source

For a spherical head with the ears at opposite points on its surface, a planar wave from a distant sound source reaches an ear on a straight line through air or if blocked by the head then partially through the shortest path along the surface of the sphere. If the sound source is at the same vertical coordinate $z$ as the ears, we can think of the head as a circle on the $x, y$ plane:

For simplicity, the head is a unit circle centered at $x = 0, y = 1$, turned so that the nose moves by an angle $\alpha$ counterclockwise from the origin. We are looking for the distances the sound still has to travel after the planar wavefront first reaches the head. There is a shortest straight line path (red dashed) of length $1-\sin(\alpha)$ to the right (red) ear. The shortest path (white dashed) to the left (white) ear consists of a straight line segment of length $1$ plus a curved path of length $\alpha$. The difference between the lengths to the two ears is $\alpha + 1 - (1 - \sin(\alpha)) = \sin(\alpha) + \alpha$. This holds for $0 \le \alpha \le \frac{\pi}{2}$, but because the function is antisymmetrical, it applies also to $-\frac{\pi}{2} \le \alpha \le 0$ giving a negative value for that range.

Modifying the equation to include the speed of sound $c$ and the head radius $r$, the inter-aural time difference (ITD) turns out as:

$$\text{ITD}(\alpha) = \frac{r}{c} (\sin(\alpha) + \alpha)$$

The formula is applicable for angles $-\frac{\pi}{2} \le \alpha \le \frac{\pi}{2}$, where $\alpha = 0$ means the source is straight ahead and $\alpha = \frac{\pi}{2}$ means it is 90° to the right, if the left channel is delayed for positive ITD. The plotted $\sin(\alpha) + \alpha$ part swings from $-\frac{\pi}{2} - 1$ to $\frac{\pi}{2} + 1$, so for a realistic maximum ITD of 0.6 ms (for distant sources) one can use simply:

$$\text{ITD}(\alpha) = \frac{0.6\ \text{ms}}{\frac{\pi}{2} + 1} (\sin(\alpha) + \alpha) \approx 0.23\ \text{ms}\ (\sin(\alpha) + \alpha)$$

For sources that are behind, before using the formula, add or subtract $\pi$ to/from $\alpha$ to bring it to $-\frac{\pi}{2}\ldots\frac{\pi}{2}$ and flip its sign.

Spherical head, nearby source

If one wants to take into account the distance to the source, the wavefront becomes spherical, or circular on the $x, y$ plane. Let's look at a head with just one ear considered:

The head is again a unit circle, but this time centered at the origin. The source (blue) is on the negative $y$ axis at a distance $d$ from the center of the head. The ear (red) is rotated by an angle $\alpha$ counterclockwise from the positive $x$ axis, or clockwise if $\alpha$ is negative. We shall constrain $-\frac{\pi}{2} \le \alpha \le \frac{\pi}{2}$ so that the shortest air path to the ear is always via the right side of the head as pictured. At the time the sound first reaches the head, the wavefront (blue) is cicular with radius $d-1$. We don't need that for the calculations, because we can use the time the sound leaves the source as the starting point. There is some negative angle $\beta$ such that if $\alpha \le \beta$, there is a direct air path of length $\sqrt{\cos(\alpha)^2 + \left(d + \sin(\alpha)\right)^2}$ to the ear located at $x = \cos(\alpha)$, $y = \sin(\alpha)$. Otherwise the sound travels a distance $D = \sqrt{\cos(\beta)^2 + \left(d + \sin(\beta)\right)^2}$ to the point $x = \cos(\beta)$, $y = \sin(\beta)$ where the straight path and the circle meet tangentially, and transitions to a further curved path of length $\alpha-\beta$. There is a right triangle (pink) that enables to solve $\beta$ using the pythagorean theorem:

$$d^2 = D^2 + 1^2\\
\Rightarrow \beta = - \text{asin}(\frac{1}{d}),\ D = \sqrt{d^2-1}$$

To recap, modifying the equations to include head radius $r$ and speed of sound $c$:

$$-\frac{\pi}{2} \le \alpha \le \frac{\pi}{2},\\
T_l = \left\{\begin{array}{ll}\frac{r}{c}\sqrt{\cos(\alpha)^2 + \left(\frac{d}{r} + \sin(\alpha)\right)^2}&\text{if } \alpha \le - \text{asin}(\frac{r}{d}),\\
\frac{r}{c}\left(\sqrt{(\frac{d}{r})^2-1} + \alpha + \text{asin}(\frac{r}{d})\right)&\text{otherwise,}\end{array}\right.\\
T_r = \left\{\begin{array}{ll}\frac{r}{c}\sqrt{\cos(\alpha)^2 + \left(\frac{d}{r} - \sin(\alpha)\right)^2}&\text{if } \alpha \geq \text{asin}(\frac{r}{d}),\\
\frac{r}{c}\left(\sqrt{(\frac{d}{r})^2-1} - \alpha + \text{asin}(\frac{r}{d})\right)&\text{otherwise,}\end{array}\right.\\
\text{ITD}(\alpha) = T_l(\alpha) - T_r(\alpha)$$

Again, for sources that are behind, before using the formula, add or subtract $\pi$ to/from $\alpha$ to bring it to $-\frac{\pi}{2}\ldots\frac{\pi}{2}$ and flip its sign.

For increasing distances, this definition of ITD approaches the previous one:


You should use py-ecg-detectors
Siply install by doing

pip install py-ecg-detectors

Then you can use for instance the well known Pan Tompkins algorithm to find the R-peaks

Here I used an ECG recording from the Apnea-ECG Database

from ecgdetectors import Detectors
import matplotlib.pyplot as plt
import pandas as pd
fs=100 # sample freq

heartbeat = pd.read_csv("/myworkdir/mitbih-database/100.csv")

detectors = Detectors(fs)

r_peaks_pan = detectors.pan_tompkins_detector(heartbeat.iloc[:,2][0:1000])
r_peaks_pan= np.asarray(r_peaks_pan)

plt.plot(heartbeat.iloc[:,2][0:1000])
plt.plot(r_peaks_pan,heartbeat.iloc[:,2][0:1000][r_peaks_pan], 'ro')

If you want you can also add a R-peak correction algorithm:

def R_correction(signal, peaks):

    num_peak=peaks.shape[0]
    peaks_corrected_list=list()
    for index in range(num_peak):
        i=peaks[index]
        cnt=i
        if cnt-1<0:
            break
        if signal[cnt]<signal[cnt-1]:
            while signal[cnt]<signal[cnt-1]:
                cnt-=1
                if cnt<0:
                    break
        elif signal[cnt]<signal[cnt+1]:
            while signal[cnt]<signal[cnt+1]:
                cnt+=1
                if cnt<0:
                    break
        peaks_corrected_list.append(cnt)
    peaks_corrected=np.asarray(peaks_corrected_list)
    return peaks_corrected

corrected_R_peak=R_correction(heartbeat.iloc[:,2][0:1000],r_peaks_pan)
plt.plot(heartbeat.iloc[:,2][0:1000])
plt.plot(corrected_R_peak,heartbeat.iloc[:,2][0:1000][corrected_R_peak], 'ro')


Question about ramp filter used in filtered backprojection

To see if your math is correct, it is useful to first understand what is, in general, the effect of a filter on a signal, and then see if you can predict what the theoretical result should look like.

If a filter has frequency response $H(Q)$, this means that its response to an input $e^{j2\pi Q_0 r}$ is the signal $H(Q_0)e^{j2\pi Q_0 r}$. In other words, a sinusoidal input of frequency $Q_0$ produces an output of the same frequency, but with amplitude $|H(Q_0)|$ and phase $\angle H(Q_0)$.

In your first question, the input has frequency $Q_0=0$. The filter's response at that frequency is $|Q_0|=0$. Then, the filter's output will be 0: frequency $Q_0=0$ is completely absorbed by the filter and it does not appear at the output.

In your second question, the input has frequency $Q_0=w/2\pi$. The filter's response at that frequency is $|Q_0|=w/2\pi$. The output, then, should be $\frac{w}{2\pi}\sin(wr)$.

As you can see, following this line of reasoning we obtain the same results as you did by just solving the equations. This corroboration should give you confidence that your results are correct, and most importantly, you can understand what the filter is doing physically.


Advise on detecting anatomical landmarks in a CT reconstructed volume

I hesitate to write this as an answer, but given that you're asking only for advice, I will do so.

I suggest investigating techniques based on the Dual-Tree Complex Wavelet Transform (DTCWT). These have shown to be useful for generating descriptors that have good tolerance to shift, scale and rotation of the source images. It's not the classic problem in that you're not allowing the points to be assigned for you, but I suspect with some thought you can adapt the techniques to predefined landmarks.

Clearly, the landmarks have some interest from the perspective of a clinician, so there is something of interest about them - it's simply a case of modelling that in the descriptor. Wavelet techniques (in particular the DTCWT) tend to be good at modelling features that the eye picks up on.

A starting point would probably be this fairly recent paper.


Equivalent of Fourier transform for Gaussians?

The (continuous) Fourier transform of a Gaussian is a Gaussian:

$$ \mathscr{F}\Big\{ x(t) \Big\} \triangleq \int\limits_{-\infty}^{\infty} x(t) \, e^{-j 2 \pi f t} \, dt $$

$$ \mathscr{F}\left\{ e^{-\pi t^2} \right\} = e^{-\pi f^2} $$

$$ \mathscr{F}\left\{ e^{-\pi \alpha t^2} \right\} = \frac{1}{\sqrt{\alpha}} e^{-\frac{\pi}{\alpha} f^2} \qquad \alpha > 0 $$

$$ \mathscr{F}\left\{ e^{-\pi \alpha (t-\tau)^2} \right\} = \frac{1}{\sqrt{\alpha}} e^{-\frac{\pi}{\alpha} f^2} e^{-j 2 \pi f \tau} $$

$$ \mathscr{F}\left\{ \sum_m c_m e^{-\pi \alpha_m (t-\tau_m)^2} \right\} = \sum_m \frac{c_m}{\sqrt{\alpha_m}} e^{-\frac{\pi}{\alpha_m} f^2} e^{-j 2 \pi f \tau_m} $$

note these gaussians have "tails" that go on forever, but gaussians usually die off to very close to zero very rapidly.

assuming no time or frequency aliasing, let $x[n] \triangleq x\big(n/f_\text{s}\big)$ where $f_\text{s}$ is the sample rate.  Then the DTFT is

$$\begin{align}
X\big(e^{j\omega}\big) &\triangleq \sum_{n=-\infty}^{\infty} x[n] \, e^{-j \omega n} \\
&= f_\text{s} \sum_{n=-\infty}^{\infty} x\big(n/f_\text{s}\big) \, e^{-j f_\text{s} \omega (n/f_\text{s})} \frac{1}{f_\text{s}} \\
& \approx f_\text{s} \int\limits_{-\infty}^{+\infty} x(t) \, e^{-j 2 \pi \frac{f_\text{s} \omega}{2 \pi} t} \, dt \\
&= f_\text{s} \ \mathscr{F}\Big\{ x(t) \Big\} \Bigg|_{f=f_\text{s}\frac{\omega}{2 \pi}} \\
\end{align}$$

so if we let $$x(t) = \sum_m c_m e^{-\pi \alpha_m (t-\tau_m)^2}$$

$$\begin{align}
x[n] &= \sum_m c_m e^{-\pi \alpha_m ((n/f_\text{s})-\tau_m)^2} \\
&= \sum_m c_m e^{-\pi (\alpha_m/f_\text{s}^2) (n-f_\text{s}\tau_m)^2} \\
&= \sum_m c_m e^{-\pi \widehat{\alpha_m} (n-\widehat{\tau_m})^2} \\
\end{align}$$

then the DTFT is

$$\begin{align}
X\big(e^{j\omega}\big) & \approx f_\text{s} \ \sum_m \frac{c_m}{\sqrt{\alpha_m}} e^{-\frac{\pi}{\alpha_m} f^2} e^{-j2\pi f\tau_m} \Bigg|_{f=f_\text{s}\frac{\omega}{2 \pi}} \\
&= f_\text{s} \ \sum_m \frac{c_m}{\sqrt{\alpha_m}} e^{-\frac{\pi}{\alpha_m} \left(f_\text{s}\frac{\omega}{2 \pi}\right)^2} e^{-j2\pi \left( f_\text{s}\frac{\omega}{2 \pi} \right)\tau_m} \\
&= f_\text{s} \ \sum_m \frac{c_m}{\sqrt{\alpha_m}} e^{-\pi\frac{f_\text{s}^2}{4 \pi^2 \alpha_m} \omega^2} e^{-j \omega (f_\text{s}\tau_m)} \\
\end{align}$$

Then, assuming you choose your DFT length $N$ to be sufficiently large (you're saying in the millions) to cover all of the gaussians with their specific delays $f_\text{s}\tau_m$ in samples and widths of ca. $\frac{f_\text{s}}{\sqrt{\alpha_m}}$.  Then the $N$-point DFT will sample this DTFT at $N$ equally spaced points around the unit circle.

$$\begin{align}
\mathcal{DFT:} \quad X[k] &= \sum_{n=0}^{N-1} x[n] \, e^{-j 2 \pi \frac{nk}{N}} \\
&= X\big(e^{j\omega}\big) \Bigg|_{\omega = 2\pi\frac{k}{N}} \quad \mathcal{:DTFT}\\
&= f_\text{s} \ \sum_m \frac{c_m}{\sqrt{\alpha_m}} e^{-\pi\frac{f_\text{s}^2}{4 \pi^2 \alpha_m} \omega^2} e^{-j \omega (f_\text{s}\tau_m)} \Bigg|_{\omega = 2\pi\frac{k}{N}}\\
&=  \ \sum_m \frac{c_m}{\sqrt{\alpha_m/f_\text{s}^2}} e^{-\pi\frac{f_\text{s}^2/\alpha_m}{4 \pi^2} \left(2\pi\frac{k}{N}\right)^2} e^{-j 2\pi\frac{k}{N} (f_\text{s}\tau_m)} \\
&= \sum_m \frac{c_m}{\sqrt{\widehat{\alpha_m}}} e^{-\frac{\pi}{\widehat{\alpha_m}} \left(\frac{k}{N}\right)^2} e^{-j 2\pi\frac{k}{N} \widehat{\tau_m}} \\
\end{align} $$

So summing up and I will remove the little hats in the discrete-time domain, if you have a collection of gaussian pulses with heights of $c_m$, peak width parameters $\frac{1}{\sqrt{\alpha_m}}$ (in units of samples), and peak positions at $\tau_m$ (also in units of samples), then your input should look like:

$$ x[n] = \sum_m c_m e^{-\pi \alpha_m (n-\tau_m)^2} $$

and the result of the $N$-point DFT should look like:

$$ X[k] \approx \sum_m \frac{c_m}{\sqrt{\alpha_m}} \ e^{-\pi/(N^2\alpha_m) k^2} \ e^{-j 2\pi(\tau_m/N) k} $$

for $-\tfrac{N}{2} < k < \tfrac{N}{2}$ .  Remember with the DFT, $X[k+N]=X[k]$ for all integer $k$.


What is Magnitude and Phase actually represents in Fast Fourier Transform?

This question is longstanding, and not fully resolved, IMO.
It hinges upon data dimensionality (1D, 2D, 3D), size (short, large), sampling (loose, tight), and mostly morphology (smooth vs. sharp).

Basically, in natural images with distinguishable objects of different contrast/texture at different scales, the data is far from stationary. Fourier applies, but sines can have very complicated relative locations, therefore phases. And to make sense of an object, sepia or luminance-distorted images still work. So, disregarding relative amplitude of components, what may matter is their relative position. Hence, phase can be primordial.

But for more stationary 2D data, with less "objects", phase can have a minor effect.


Ok so here is what I found. The distance is dependent on the way that the mfccs are calculated. This makes sense to me and also explains why cepstral mean normalization affects the values of the MCD.

I found this implementation (https://github.com/MattShannon/mcd), which unfortunately did not support .wav files. I ran this and it gives results that are in the correct range for its test examples. I looked into the code and found that the the difference was not so much in the calculation of the MCD (which is straight forward) but rather the mfccs. The library I was using gave mfccs at a different scale (I tried python_speech_features and pyspeech). Using the mfccs from Matlab gave me mfccs in a similar range to that of the repo so I opted to use those.

A noteworthy point is that the first coefficient is quite large which is why it is ignored in (https://github.com/MattShannon/mcd) so I did the same. There are some artifacts in the speech I synthesise which likely accounts for why my MCD is still above 10 but its still lower than 30 now.


This is probably the craziest, though important question on DSP SE ;) Answer is, because they:

look like peak:

and tend to have a shelf either in the beginning or the end of frequency range:


GnuRadio: Clock Recovery MM: imu out of bounds

The M&M clock recovery block has a bug in it, such that it can't reliably handle more than 8 samples per symbol without inadvertently dropping samples.  That interpolating resampler alignment bug may also be responsible for generating mu values outside of [0.0, 1.0].

I'd strongly recommend, downsampling the input from 2.5ksamples/symbol down to 6 or 8 samples /symbol before the M&M block.


What would be the proper notation for denoting part of signal that satisfies a condition?

What you are looking for is called an Indicator Function which would assign a $1$ for all those elements of $E(t)$ that satisfy the given condition.

The same could also be expressed via the use of Set Builder Notation as $E_{interested} = \{1|e \in E \land e>C\}$ translating to "Elements of $E_{interested}$ take the value of 1 for all elements $e$ in $E$ satisfying the condition $e>C$". You might have to add $E_{interested} = \{0|e \in E\}$ as a pre-initialisation step to show the default value of $E_{interested}$.

The difference between the two is that the indicator function is defined and so you can simply reference it (e.g. "...where $\mathcal{\chi}(E, \Phi(E))$ denotes the Indicator Function using a predicate $\Phi$...", or even drop $\Phi$ if the predicate is the same throughout) while the set builder notation version could be used more as a "one-of" in-line way in writing down an algorithm.

To an extent, it is also a matter of notation already used and what you are trying to convey.

Hope this helps.

EDIT :

Just to clarify, if you are also interested in the actual values of $E(t)$ (interrupted by zeros wherever $E(t) \leq C$), then given $E_{interested}$ you can now $E_{interested} \cdot E$ to produce that sequence that contains the actual values of $E(t)$ satisfying the condition (rather than just the "pulse sequence").


Camera calibration using multiple images

You are correct: to calibrate a camera you need a correspondence between 3D world points and 2D image points. The problem is that the 3D points cannot be co-planar, so people were building 3D calibration rigs, e. g. a box made of checkerboards. One image of a rig like that would be enough to calibrate, but those rigs are hard to build, because you have to get the planes to be at exactly right angles to each other.

Then Zhengyou Zhang came up with the algorithm to calibrate from multiple images of a planar pattern, like a checkeboard. You still need non-coplanar points, but you get them from multiple images of a planar pattern in different 3D orientations, rather than a single image of a 3D rig.

Noise and errors are another reason why you need multiple images. You cannot detect the checkerboard corners with infinite precision. So you need a lot of points, to estimate the camera parameters more reliably.


Processing Accelerometer Vibration Data With FFT

After taking a closer look at your results, here is what I think.

Because you are playing back sinusoid through your speaker you are also observing it's properties. It has some significant harmonic distortions appearing as the harmonics of fundamental frequency: $200$, $300$, $400 \mathrm{Hz}$ in your case. That is either because:

You are driving your speaker at high volume
It is simply crap for accurate measurements (sorry...)

Obviously you could also get harmonic distortions coming from your sound card, but this is unlikely.

On the other hand what you see at $50 \mathrm{Hz}$ is the mains hum of AC voltage (you are definitely not in USA). This distortion also has it's harmonics at $100$, $150$, $200 \mathrm{Hz}$, etc. So you can see that they do overlap with your test signal.

Nevertheless if you look at higher frequencies, i.e. around $400 \mathrm{Hz}$, you will see two peaks next to each other. Those are two harmonics, one of your test signal and other one of your mains hum. Because there was a slight offset at low frequencies it diverges even more at higher frequencies.

In order to confirm the mains hum you should simply record silence. If it is present, then well - you have a slight problem. If you are using laptop, try to unplug the AC adaptor - it should help. Otherwise you are simply picking up everything with your power supply. You can also try external audio interface. On the contrary, if you are already using some kind of interface with XLR microphone or whatever - try to break the ground loop.

When it comes to speaker distortions try to use other one or maybe use it a slightly lower volume.


How to convert dB back to manitude and then to rectangular format

If $X_\mathrm{db} = 20 \cdot \mathrm{log}_{10}(|X|)$, then $|X| = 10^{X_\mathrm{db} / 20}$, of course. Furthermore, you can compute the cartesian representation of a complex number $c = a + ib$ from its magnitude $|c|$ and phase $\varphi$:  $c = |c| \cdot \mathrm{exp}(i \varphi)$

So you are probably looking for something like this:

import numpy as np
from matplotlib import pyplot as plt

fs = 100
t = np.arange(100)
x = np.sin(2 * np.pi * 3 * t / fs)

X = np.fft.fft(x)
X_mag = np.abs(X)
X_db = 20 * np.log10(X_mag)
X_phase = np.angle(X)

# ... do modification on X_db and X_phase ...

X_mag_rec = 10 ** (X_db / 20)
X_rec = X_mag_rec * np.exp(1j * X_phase)
x_rec = np.real(np.fft.ifft(X_rec))

plt.plot(t, x, label='original')
plt.plot(t, x_rec, label='reconstructed')
plt.legend()
plt.show()

If you are doing the modification in the STFT domain, note there may be no valid inverse STFT. For details on that look into:

D. Griffin, J. Lim: Signal estimation from modified short-time Fourier transform, in: IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 32, no. 2, pp. 236–243, Apr 1984.


When applying a bandpass filter designed in Matlab's DSP system toolbox, is it necessary to do any normalization?

Note that you are cutting a lot of signal frequencies in the transition from 10Hz to 200Hz (look the attenuation in the frequency response of the filter), consider to use lowpass filter or to change filter prototype to someone with faster transition.
By cut definition each frequency that is attenuated more than -3dB is cut.


3D Shape Nets is published this year at CVPR (2015). You could access it here. This is a good way to use voxels as input to DNNs. There are also works applying DNNs to medical imaging such as this one and this one. I guess these would help you a lot.


Is multiplying by Ts faster than dividing by fs?

Generally, it makes sense to ensure that your code is logically correct, that it is numerically well-behaved, intuitive to read and tested. That is hard enough. Only when you observe that some innerloop or library call is a real hotspot, affecting the functionality of your software does it make sense to rewrite code for speed, and then you should always profile before and after.
If a constant is known compile time, the compiler may apply the inversion to substitute division for multiplication, if this is within precision constraints and runs faster for a given target. If ppssible, I would rather outsource that complexity to the compiler.
Edit:
It is not "technically the same operation". To see why, have a look at this MATLAB snippet:
a = single(10.0)
b = 1/a
c = 42/a
d = 42*b

c-d

ans =

  single

 -4.7684e-07

Since floating-point is operating with finite precision and intermediate rounding, the order of operations does matter. Depending on compiler flags, the compiler may be allowed to re-order floatingpoint arithmetic even though the result will differ to some degree.
If we look at the binary representation we see that they differ in the lsb:
dec2bin(typecast(c,'uint32'),32)
dec2bin(typecast(d,'uint32'),32)

ans = '01000000100001100110011001100110'
ans = '01000000100001100110011001100111'

-k


Detecting a low frequency cue in audio recording

You can solve this problem with a simple time based low pass filter, no need for an FFT:

Calculate the minimum integration time you want by dividing your sampling rate by the frequency you want to detect. For a 10Hz at 44100Hz, it would be 4410 samples, so let's take 4500 samples.
Then apply a low pass filter to this group of samples (choosing the right frequency cut), you will isolate the frequencies you want. I do not really know how to use them, but it should look something like this:

$$y(n) = x{(n)}+x(n-1) + ...$$

Calculate the relative power of your filtrated group of samples in dBFS (Full Scale) based on this formula:

$$p_{RMS} = \sqrt{\frac{ x_1^2 + x_2^2 + \ldots}n } $$
$$dbFS = 20\log_{10}\frac{ p_{RMS}}{p_{max}}$$
'n' behind the number of samples and 'pmax' the maximum value of a sample.

If the level you get is superior to a specific threshold of your choice (not to low, because your audio certainly have a tiny low frequency component), then there is a cue, otherwise there isn't.

As for the timestamp aspect, you can increase a variable each time you get into this low frequency detection loop, and calculate the time passed using the sampling rate. Then write an array containing the input/output time of each successful detection.

Not sure if the minimum integration time is really necessary though.


Here is a simple way of looking at it which avoids any algebra and uses minimal symbolism. Let us calculate the output $y$ for an input $x$ which is a unit pulse function. We have that

$y[-2] = (-2)(x[-2])^2 = 0$ since $x[-2] = 0$.
$y[-1] = (-1)(x[-1])^2 = 0$ since $x[-1] = 0$.
$y[0] = (0)(x[0])^2 = 0$ (regardless of the value of $x[0]$).
$y[1] = (1)(x[1])^2 = 0$ since $x[1] = 0$.
$y[2] = (2)(x[2])^2 = 0$ since $x[2] = 0$.

and so on.  The entire input is blocked and the output $y$ is always $0$.

What about the output $\hat{y}$ when the input $\hat{x}$ is just the previous $x$
delayed by one time unit? The only difference now is that $\hat{x}[1] = 1$ and
we have that

$\hat{y}[1] = (1)(\hat{x}[1])^2 = 1$ since $\hat{x}[1] = 1$.

It is, I hope, obvious that $\hat{x}[n] = 0$ for all integers $n$ except $n=1$.
So, since the input $\hat{x}$ is just a delayed version of input $x$, is
is true that the output $\hat{y}$ is just a delayed version of output $y$?
What, if anything, does this tell you about the system properties?


Is a change in colorspace needed when comparing colors only evaluated by the computer?

Delta E color difference formula is based on CIELAB color space and human perception. To convert from RGB to CIELAB space, you would need a base whitepoint (standard illuminant) which in this situation you do not have and need.

Your assumption is correct. You do not specifically need to have a color difference similar to human perception. Additionally color distance metric should be based on your needs (you can use weighting depending on your needs).


What is the correct way to apply multiple IIR filters to a signal in MATLAB

If you filter out all frequencies above a certain cut-off frequency by a low pass filter, and if the filter is good, then there's hardly anything left above that cut-off frequency. If you then filter out everything below a certain cut-off frequency using a high pass filter, and if that cut-off frequency is greater than the cut-off of the low pass, then you eliminate the rest. There's not much left to do for the band pass filter afterwards, because you've basically suppressed your whole signal. That's why you don't hear anything.

What you probably wanted to do was apply the filters separately, i.e. create one low pass output, one high pass output, and one band pass output. In this case you have to filter the original signal x with the three different filters to generate three different output signals.

yl = filter(lowpass,x);
yh = filter(highpass,x);
yb = filter(bandpass,x);


Fourier transform possible on non-rectangular part of an image

Even when the image part is rectangular, the FFT result will be contaminated by boundary effects. The FFT of the part is the Fourier transform of an infinite tiling of the rectangular part. (An elliptical or circular part cannot be tiled. A hexagonal tiling would get pretty close though.) If there is a discontinuity going from one tile to another, this will create high-frequency components that are not descriptive of the insides of the part. Randomizing the phases will spatially translate these components so that they appear also inside the part, which is probably not what you want. Mirroring at the boundaries does not eliminate discontinuities of the derivative, and requires doubling of the FFT dimensions.

This is where a smoothly tapering 2-d window function can help. For a circular part you might like the Tukey window with the distance from the part center as its argument.

For an ellipse, stretch the window function accordingly. Multiply the data by the window function to select the area you want to analyze, and FFT the result. Then do the phase randomization. Mask the noise by multiplying it with an ellipitical disc-like mask that has value 0 outside the disc and value 1 inside the disc. You can also use the window function as the mask to get a fuzzy boundary. I did not find yet a perfect way to normalize the noise, but matching the mean and the standard deviation both weighted by the mask works quite OK:

Multiply the input image by the complementary mask: 1 - mask, and sum the two together. Octave source:

pkg load signal;
r = 120; # Circle radius
n = 512; # Image width and height
tukeyr = 0.5; # Window transition relative size
lena = double(imread("lena.gif"))/255;
tukeyn = round(r/(1-tukeyr));
x = repmat([1:n], n, 1);
y = repmat([1:n]', 1, n);
dist = sqrt((x-n/2).^2 + (y-n/2).^2);
tukey = tukeywin(tukeyn*2, tukeyr)(tukeyn+1:tukeyn*2);
win = reshape(interp1(1:tukeyn, tukey, reshape(dist, n*n, 1), 0), n, n);
randphase = 2*pi*rand(n, n);
noisefft = abs(fft2(lena.*win)).*exp(j*randphase);
noise = real(ifft2(noisefft));

#Circular mask
weight = r-dist;
weight = (weight >= 1)*1 + (weight > 0 & weight < 1).*weight;
#Or use window as mask
#weight = win;

weightsum = sum(sum(weight));
noiseweightedmean = sum(sum(noise.*weight))/weightsum;
noiseweightedstdev = sqrt(sum(sum((noise-noiseweightedmean).^2.*weight))/weightsum);
lenaweightedmean = sum(sum(lena.*weight))/weightsum;
lenaweightedstdev = sqrt(sum(sum((lena-lenaweightedmean).^2.*weight))/weightsum);
normalizednoise = ((noise - noiseweightedmean) / noiseweightedstdev * lenaweightedstdev) + lenaweightedmean;
imshow(lena.*(1-weight).+normalizednoise.*weight);

P.S. I found that highpass filtering the copy of the Lena image used for noise generation, by setting some of its lowest-frequency bins to zero, will reduce the error in the magnitudes of the fft of the final image compared to fft of the source image. I think the reason is that low-frequency amplitude fluctuations will be reduced. Within the circle (ellipse) those appear in practice as zero frequency, and in their absence less normalization is needed. So that's one thing to investigate.


You can get the sound pressure level from PSD, and the velocity is related to the sound pressure. According to the equation of motion,
$$
\rho \frac{\mathrm{d}\vec{v}}{\mathrm{d}t} = -\nabla p \tag{1}
$$
where $\rho$ is the density of the medium, $\vec{v}$ is the particle velocity, $\nabla=\frac{\partial}{\partial x}\vec{i}+\frac{\partial}{\partial y}\vec{j}+\frac{\partial}{\partial z}\vec{k}$ is the gradient operator, and $p$ is the sound pressure. When the sound is not so loud, the equation of motion can be linearized and the total derivative $\mathrm{d}\vec{v}/\mathrm{d}t$ becomes partial derivative $\partial{\vec{v}}/\partial{t}$.
Furthermore, if you make a plane wave assumption, then the particle velocity amplitude is
$$
v = \pm \frac{p}{\rho_0 c_0} \tag{2}
$$
for plane waves of forward and backward propagation, respectively. $\rho_0$ and $c_0$ are respective the medium density and sound speed without acoustic disturbance. Their product $\rho_0 c_0$ is the characteristic specific acoustic impedance of the medium.
However, your data is acquired in a cavity with an extremely high level of PSD, thus the linearization is obviously inappropriate. The total derivative in Eq. (1) should be
$$
\frac{\mathrm{d}\vec{v}}{\mathrm{d}t} = \frac{\partial \vec{v}}{\partial t} + \frac{\partial\vec{v}}{\partial x} \frac{\partial x}{\partial t}+ \frac{\partial\vec{v}}{\partial y} \frac{\partial y}{\partial t}+ \frac{\partial\vec{v}}{\partial z} \frac{\partial z}{\partial t} \\
= \frac{\partial \vec{v}}{\partial t} + v_x\frac{\partial\vec{v}}{\partial x}+ v_y\frac{\partial\vec{v}}{\partial y}+ v_z\frac{\partial\vec{v}}{\partial z}
\tag{3}
$$
In addition, velocity is also related to frequency, but it seems that you want an overall velocity. You can obtain the time-domain sound pressure according to the Parseval's theorem, and then calculate the particle velocity.


Importance of Discrete Cosine Transform over High Pass filter?

Yes, if you apply a filter, say a high Pass, you get the low frequency components of the image, but still in the image domain. You need a Transformation, like the DCT or the FFT to get to the frequency domain, to see the separate frequency components of the image.


RMS from calculation from PSD using Welch's method

Also, you can use the PSD for estimation of the RMS, the easier way for estimation of the RMS is a direct calculation.

In case, you still want to use the PSD for estimation of the RMS, you should normalize the values by the square of the average value of the window. For the Hann window, the normalization factor will be $\frac{1}{4}$.

update
The question is a duplicate from RMS calculation in frequency domain after windowing


How to detect lines arranged in a visual pattern in an image?

If it is a true bw-image, i.e. it contains binary values only, you can use the Hough transform for circles (or a special one for ellipses). In the Hough space you will find maxima for the coordinates of the center and the radius or radii, respectively. If necessary, you can use a statistical test to proof the maxima since there are offen multiple local maxima.


Number of datapoints after decimation in discrete wavelet wavelet analysis

As with Fourier, as with IIR filtering, discrete wavelet forward implementation, and how they can be inversed, depends of lot on how you treat the samples that you don't know, for standard finite support data. In other words: as you convolve, what do you do on the left side and the right side of the signal.

There is a huge literature on how to perform a non-expansive decomposition, that generates the same number of coefficients that the number of samples. But yo can see the constraints easily: when I have 3 samples, how to subsample this evenly into a low-pass and a high-pass series?

To made a long story short, if the signal length can be divided by $2^L$, $L$ wavelet levels can be obtained in general with the appropriate signal extension, governed by Matlab's dwtmode: "the options The DWT associated with the symmetric, smooth, zero, and periodic extension modes are slightly redundant."
This is generally due to how wavelet filters overlap where the signal is not defined. This is related to the convolution with lowpass and highpass filters as well, before decimation. And convolution tends to increase the signal's length. A filter with support $K$ convolved with a signal of length $N$ generally yields a filtered signal of length $N+K-1$. So if you decimate and when to keep all necessary information, you need to keep a little more coefficients at each level. There are nice description in Chapter 4 of "Wavelets in chemistry", Beata Walczak, 2000.

However in some case, like with the Haar wavelet, or as detailed in C. Brislawn, 1996 Classification of Nonexpansive Symmetric Extension Transforms for Multirate Filter Banks, nonexpansive schemes can be designed.

In practice, people who really want to be nonexpansive uses modified DWT algorithms, like the lifting scheme, or wavelets on the interval, that are not 100% pure discrete wavelet transforms.


You can use partial fractions first to expand the $X(e^{j\omega})$ as

$$X(e^{j\omega})=\frac{\tfrac{7}{4}}{1+\frac{1}{2}e^{-j\omega}}+\frac{\tfrac{5}{4}}{1-\frac{1}{2}e^{-j\omega}} $$.

Then taking the inverse DTFT you get $\tfrac{7}{4}(-\frac{1}{2})^{n}u(n)+\tfrac{5}{4}(\frac{1}{2})^{n}u(n)$


Bilinear is the most widely used method. The nearest neighbor down-sampling algorithm is the fastest but least accurate.

Note that when you are trying to down-scale an image by half, the bilinear sampling algorithm becomes the average (mean as you mentioned):

x1[i/2][j/2] = (x[i][j] + x[i+1][j] + x[i][j+1] + x[i+1][j+1]) / 4;


Mathematical advantages of the ZT, DTFT and DT?

My friend, i hope you can gain a proper command of those techniques for achieving an excelent result at your exam. Of course the question is broad and more discussionfull, but the purpose of this is "clarify" not to turn this into a "closary" of questions.

Indeed. there are a lot of variants for the transforms, but only a few of them have a clear direct meaning, and the others are generalizations seldom used.

At first, the standard, classical FT, is the main entry for your analysis.
$$\mathbf{F}=X(f)=\int_{-\infty}^{\infty}x(t)e^{-i2 \pi ft}dt$$
$$\mathbf{F}^{-1}=x(t)=\int_{-\infty}^{\infty}X(f)e^{i2 \pi ft}df$$
That is the base, plain and simple, sweet, the FT and the IFT are the same (of course exchanging the signs and variables :) ) and with the easiest expressions $\mathbf{F}\{e^{-\pi t^2}\}=e^{-\pi f^2}$, $\mathbf{F}\{1\}=\delta(f)$, $\mathbf{F}\{\cos(\pi t)\}=\frac{1}{2}\delta(f-1)+\frac{1}{2}\delta(f+1)$, $\mathbf{F}\{\Delta(t)\}=\Pi(f)$, $\mathbf{F}\{\Pi(t)\}=\Delta(f)$, $\mathbf{F}\{|||(t)\}=|||(f)$ (where $\Delta$, $\Pi$ and $|||$ are the triangle, step and comb functions:)), and so on.
You use this on the handwritting, papers, reports, analysis, etc.

We can safely ommit the DTFT and DFFT, and the tricky $i$ and $k$ symbols for discrete time and frequency and keep $t$ and $f$, discretized, with the sampling time $dt$  and the frequency periodic spacing $df$ (again with the same notation as the integral). Thus for the DTFT direct transform is:
$$\mathbf{F}_{DT}=\sum_{t=-\infty}^{\infty}x(t)e^{-i2\pi ft}dt$$
and the DFFT inverse is:
$$\mathbf{F}^{-1}_{DF}=\sum_{f=-\infty}^{\infty}X(f)e^{i2\pi ft}df$$

because the $|||(t)$ (and $|||(f)$) comb function:
$$|||(t)=\sum_{i=-\infty}^{\infty}\delta(t-i)$$

is all what we need for recovering the related expressions and cumbersome hard to learn theorems regarding sampling (DT) and perioditizicing (DF). And you recover all your symmetries just by multiplying by your combs:

$$\mathbf{F}_{DT}(x(t))=\mathbf{F}(x(t)\frac{1}{dt}|||(\frac{t}{dt}))=X(f)\frac{1}{dt}\frac{1}{F}|||(\frac{f}{F}), F=\frac{1}{dt}$$

Hence, at my second place, you have your well forgotten ZT, unilateral, again considering the $t$ as a discretized variable under a sampling time $dt$:
$$X(z)=\sum_{t=0}^\infty x(t)z^{-t}$$

You have a DAQ with a lot of sensors and controllers, and you have some fixed sampling time and you need to know the continuous parameters of the system based on the discrete behaviour, or tuning some controllers, or identifying a system. Soon or later you need to pass onto the Z domain and check if your IIR filters are stable or nor, or if your FIR smoothers are too short on delay, and comparing your ARMA models with your differential equations.


How to represent the following system function in a block diagram?

It is a simple first-order recursive system like this:

In your case $b[0]=1$, $b[1]=0$, and $a[1]=-0.3(1+j)$. Of course the coefficient is complex if your system is complex. What you can do is split the real and imaginary parts and get two real-valued systems. The output of one is the real part of the output signal $y[n]$, and the output of the other is the imaginary part:

$$y_R[n]=0.3y_R[n-1]-0.3y_I[n-1]+x_R[n]\\
y_I[n]=0.3y_R[n-1]+0.3y_I[n-1]+x_I[n]$$

where $y_R[n]$ and $y_I[n]$ are the real and imaginary parts of $y[n]$, respectively (and the same for $x[n]$). If the input signal is real-valued then of course $x_I[n]=0$.

As for general complex-valued systems also have a look at this figure from the textbook Digital Communication by Lee and Messerschmitt:


How to obtain mean Power Spectrum from two audio channels?

The sound pressure level (and derived metrics such as the SPL spectrum) are dependent on the position of the listener, since it is a measure of the pressure of the wavefront moving through space. As such it doesn't make all that much sense to try to combine SPL's from different positions into a single value, except perhaps if the listeners were exactly the same distance from a single omnidirectional audio source and there was no possibility of reflections.

The equations you have posted are for two audio sources, not two listeners, so I am not sure they are relevant.

The only value of any use I can imagine you being able to calculate from two listeners would be the peak SPL / spectrum, this would simply be the maximum value of both channels for any given sample.


In general there is no systematic way and you simply have to analyze the given system. In the case of the system in your question, it's easy to see that it can't be invertible, because the output is just a constant, namely the integral over the input function (assuming this integral exists). There are infinitely many functions which will result in the same value when integrated from $-\infty$ to $\infty$.


Finding the system output by convolution

You can't compute the output by convolving the input with the system's impulse response if the initial conditions are non-zero. The reason is simple: there is no information about the initial conditions in the impulse response or, equivalently, in the system's transfer function.

What you can do is use the following properties of the unilateral Laplace transform:

$$\mathcal{L}\{y'(t)\}=sY(s)-y(0^+)\tag{1}$$

$$\mathcal{L}\{y''(t)\}=s^2Y(s)-sy(0^+)-y'(0^+)\tag{2}$$

If you use $(1)$ and $(2)$ when you transform the differential equation, you'll get an equation for the Laplace transform of $y(t)$ which takes into account the given initial conditions. The inverse Laplace transform of that expression gives you the desired system output for non-zero initial conditions.

As for stability, causal systems with a rational transfer function are stable if all poles (i.e., the zeros of the transfer function's denominator) are inside the left half-plane.


How do I apply a binary mask and STFT to produce an audio file?

The idea of the short-time Fourier transform (STFT) in this case is to compute a representation of the input signal where we see how the frequency content of the signal evolves over time. If you aren't familiar with the STFT, I would suggest you read up on it a bit first. The Wikipedia article on the subject should get you started.

To your question. The STFT function in your code gives you a complex-valued matrix in the output. When you move up or down a column in the matrix you are moving in the frequency direction. When you move along a row in the matrix you are moving along the time direction. The magnitude of a single complex-valued element in the matrix gives you an estimate of the energy in a certain frequency band at a certain time position in the input signal.

I recommend you try to understand the parameters that are given to you in the function call. You should try to understand how the window size, hop size, number of FFT bins, and the windowing function used affect the output of the function.

If we denote your complex valued STFT with $$X = \textrm{STFT\{x(n)\}},$$ and the binary mask with $M$, the output signal in the time-frequency domain is given by:

$$Y = X \odot M,$$
where $\odot$ refers to element-wise multiplication. Then you want to convert your modified time-frequency representation $Y$ back into to a time signal. This is done by doing the inverse STFT:

$$y(n) = \textrm{ISTFT}\{Y\}.$$

This should get you going, but time-frequency processing is a complex topic and I highly recommend trying to find resources to help you really understand what is going on when you are doing things such as applying a binary mask in the STFT-domain.


MSER (Maximally stable extremal regions) are regions, not points. And they're invariant to affine transformation. But it's not a segmentation method, strictly speaking

Informally speaking, the idea is to find blobs at various thresholds, then select the blobs that have the least change in shape/area over a range of thresholds. These regions should be stable for a large range of grayscale and geometric transformations.


How can I insert message bits into the frequency spectrum of an audio signal using Matlab?

The FFT is just an implementation of the DFT. There's, mathematically, no difference.

The DFT is a linear operation. So it makes absolutely zero difference whether you add up two signals in time or frequency domain.


Auditory processing has a latency of around 40ms from stimulation at the cochlea to awareness in the brain. Those 40ms appear in a number of auditory cognitive effects and mark the separation between the perceived "now" and the past. Within these 40 ms the brain can group information and produce masking effects. Pre-masking is much shorter than post-masking, because of causality, but it's not necessarily 0.


The SINR characteristics of a DSSS Multiuser Communication are very difficult to approximate theoretically.

For the Downlink at least, it is easy enough to characterize it with lesser parameters, but still difficult to obtain any particular derivation. It depends heavily upon the Spreading Sequences(their auto-correlations and cross-correlations), symbols of the interfering users(their average Energy) and so on. As a very decent approximation, $SINR$ is usually computed using random spreading sequences assumption. According to Andreas Goldsmith's textbook "Wireless Communications", for a $K$-user system with $N$ chips per symbol, the $SINR$(assuming Random signal transmissions for every user, again a probabilistic assumption) is given by $SINR = (\frac{N_{o}}{E_{s}}+\frac{K-1}{G})^{-1}$ where $G=N$ approximately.
For the Uplink it is even more a difficult task to analyze this. If all users' signals are received with the same power($AWGN$), it can be shown that for those Asynchronous users assuming Random Spreading codes with $N$ Chips per symbol, random start times and random carrier phases, $SINR = (\frac{K-1}{3N}+\frac{N_{0}}{E_{s}})^{-1}$.

For simple simulations, usually people assume that the $SINR$ can be calculated as, $SINR = \frac{E_{s}}{N_{0}+\sum_{i=1}^{K-1}\frac{E_{si}}{G}}$ where G is the processing gain which can be roughly taken to be as $K$. This factor of $G$ applies if the codes are perfectly orthogonal and its approximation to $N$ is usually justified assuming the properties of Spreading codes.

In all the above formulae, $N_{0}$ is the Noise Power which is equal to $2*\sigma^{2}$
for 2-Dimensional Signalling and $E_{s}$ is the average Received Energy of the signal. You can refer to Goldsmith if you to learn any further on the math behind it all. Now you can straight away plugin your numerical value on the LHS of the expressions that are stated above and you can go about calculating the Thermal Noise Power since you know the Received Powers $E_{s}$ for every user's signal. And all these formulae are for $AWGN$ conditions. The difficulty comes from the fact that Interferers cannot be approximated as Random Gaussian Noise.


general n-dimensional name of Pixel and Voxel

From Wiktionary:

  hypervoxel (plural hypervoxels)

  A multidimensional analogue of a voxel

From Digital Design Media - Mitchell & McCullough, 1995:


calculating the stop-band edge frequency of a low-pass FIR

The question is what you really want to do. Because the answer to your question is that there is no way to exactly compute the stopband edge from the three values you have. Looking at your specs I guess you're talking about a Parks-McClellan equiripple FIR filter, for which there are heuristic formulas relating the width of the transition band and the maximum ripple to the filter order (see this answer).

Usually in a design problem the passband and stopband edges would be given, and the maximally allowable error (ripple) would be given also, and you would need to figure out (by trial and error and/or by some heuristic formulas) what the required filter order would need to be. The problem you asked about doesn't occur in practice.

As a last comment, if you were given the actual filter coefficients then, by definition, the stopband edge is of course the lowest frequency for which the filter's magnitude response equals the given stopband ripple value (for a low pass filter, that is).


How do i select my frequency range to plot when fft is made on time domain

The Discrete Fourier Transform takes a signal from the discrete time domain to a discrete frequency one. In the process, it is possible to evaluate the DFT in less (or more) "frequency bins".

In the time domain, the signal is indexed by a variable that denotes discrete samples. For example $x[n]$ where $n$ relates to physical time by $n \cdot T_s$ where $T_s$ is the sampling period.

In the frequency domain, the spectrum of the signal is indexed by a variable that denotes discrete frequency. For example $X[k]$ where $k$ relates to physical frequency by $k \cdot \frac{Fs}{N_F}$ where $N_F$ is the length of the DFT.

So, given a signal that was acquired at $Fs = 44100$ Hz in a buffer that was 1024 samples long (that is the $x[n]$) which was then transformed by a 512 point DFT (that is the $X[k]$), the correspondence between $k \in [0 .. 511]$ and physical frequency $f$ would be $f_k = k \cdot \frac{44100}{512}$ or a spacing of 86 Hz.

You can then use this correspondence to choose which part of the spectrum you want to plot.

As a side note, your accelerometer board is pausing data logging to write the acquired buffer to an SD card or send it over some interface elsewhere. This creates very short zero value gaps, which, overall, are going to distort your spectrum into something that looks "wavy" in a very distinct way. In fact, your spectrum is effectively the spectrum of your vibration data as modulated at 100% by a square wave at whatever frequency is implied by the length of the buffer and the speed of the SD Card / interface of your data logger. Depending on how short your recordings are going to be and how they are timed, you are going to have to do something about this, otherwise your spectrum will not be the "real" one.

Hope this helps.


Why is the sound field intensity due to $K$ point sources given by $ I(p,\omega) = \sum_{k=1}^K \sigma_k^2(\omega) \delta(p - p_k)$?

I believe the confusion comes from the fact that equation $(*)$ is referring to a received intensity (in units of W/m^2), whereas the second equation is referring to a radiated intensity.

The assumption behind equation $(*)$ is that the received power as a function of angle is the surface integral of the incident intensity over some aperture area $A \subset S$. Because the surface area element $dA$ has units of area, the result of the integral will be in units of power (W).

The Dirac delta distribution $\delta(p - p_k)$ expresses that the power density is coming from exactly one direction $p_k$ and nowhere else. In this context, the Dirac delta distribution possesses the following property:
$$
\int_A \delta(p - p_k) dA' = \begin{cases} 0 & p_k \notin A \\ 1 & p_k \in A \end{cases},
$$
which is analogous to a "density" where all the "mass" is concentrated at one point $p_k$.

To calculate the total received power $P$, you integrate the received intensity over the receive aperture surface $A$:
$$
P = \int_A I(p,\omega) dA'
$$

Each term of the integral will look something like this:
$$
\int_A \sigma_k^2(\omega) \delta(p - p_k) dA' = \begin{cases} 0 & p_k \notin A \\ \sigma_k^2(\omega) & p_k \in A \end{cases},
$$
so in other words, you are just adding up the power coming from each source that your receive aperture $A$ captures, and each source contributes power $\sigma^2_k(\omega)$.

Implicit in this formulation is that all the sources are incoherent, i.e., the waves do not have a fixed phase relationship and thus the average power of the superposition is the sum of the contributions from the individual sources.


What is the difference between convolution filter, low-pass filter, and median filter?

A median filter is most certainly not a "blur" filter, purely on the basis that it tends to preserve edges. Edges are abrupt transitions of brightness and therefore that information is encoded in the high frequencies of the spectrum. Incidentally those high frequencies are the ones that low-pass filters suppress the most, leading to that "blurry" appearance because only the form of the depicted objects is retained rather than their details.

Convolution can only be used to represent linear time invariant systems. That is, systems whose relationship between the input and the output is a linear combination and therefore proportional inputs produce proportional outputs. Furthermore, this relationship does not change with respect to time.

A median operator works by sorting the values of the $M \times N$ mask of pixels surrounding some $i,j$ pixel and then assigning their median (the pixel value that happens to lie at the midpoint of the range of values) to the $i,j$ pixel.

This operation, of sorting and then selecting the median, cannot be expressed as a linear combination of the pixel values within the $M \times N$ mask, hence the name "non-linear filter".

Hope this helps.


Types of rounding in coefficients quantization

The various rounding methods have a computation vs. quantization error tradeoff.

Truncate

Truncation is the simplest method.  Everything after the decimal point is simply lopped off.  For instance, both 2.1 and 2.9 become 2.  This is very simple, but is the worst method in terms of quantization error.  It is particularly bad because when you are dealing with non-negative numbers it introduces a strong negative bias.  In some algorithms that can be bad.

Round/Convergent Rounding

Simply saying "Round" doesn't tell you enough to know what is meant.  What kind of rounding?  People usually mean convergent rounding when they say "round", so I will assume that that is what is meant.

Convergent rounding rounds down when the decimal place is $\le .499\overline{9}$ and rounds up when the decimal place is $\ge .500\overline{0}1$.  The question remains- what to do when you have exactly $.5$?  Some rounding algorithms always round down (this is called "round-to-zero"), but that introduces a very small amount of bias.  Convergent rounding tries to eliminate the bias by rounding to the nearest even number, the assumption being that around half the time that will be up and half the time it will be down.

This is the best algorithm in terms of quantization noise, but is the most computationally intense.  In many situations though, you don't care how computationally intense it is since the calculation isn't done at run-time.

Round-to-Zero

As mentioned previously, round-to-zero is the same as convergent rounding except it always rounds towards zero when the decimal portion is $.5$.  This is a compromise between truncation (easy computationally) and convergent rounding (low quantization error).


Building a FIR filter from an arbitrary frequency-magnitude response curve (eg. Least Squares fitting)

Am I understanding this correctly, and if so, will that give me what I
  need?

Yes you are understanding everything correctly. If your "need" is a filter that will match the magnitude response then all choices should work. If you are interested in a very quick and move-on solution, then proceed with Frequency Sampling with a very large number of taps and move on. If you actually have a concern with efficiency (meaning you will be implementing the filter in hardware), then I recommend Least-Squares (although you could run into convergence issues). If your interest is education, then try all and compare the results yourself as nothing will beat the first-hand experience. All approaches will result in a solution.

The descriptions of the different approaches for designing FIR filters as stated by the OP are reasonably accurate, and given no constraint on total number of taps (and total delay) all could achieve the desired objective, so often the decision is just a trade of how much time you want to spend on designing the filter to achieve the best efficiency (least number of taps).

A fourth approach not mentioned is the minimax solution (equiripple) using the Remez exchange algorithm is also a common design approach.

The graphic below lays out my flow-chart for the typical steps in the FIR design process, with further details specific to estimating the order of the filter and the 4 major design approaches.

Estimating Order: Estimating the number of taps required in the FIR filter is determined before any further design efforts are made in any of the approaches. The number of taps required is driven from the derivative of the filter response versus frequency, both in magnitude and phase. The derivative of the magnitude response for typical low pass or high pass filter structures translates to very tight transition bands require a large number of taps, and fred harris' rule which is my favorite due to it's simplicity:

$$N \approx \frac{A}{22}\frac{2\pi}{\Delta\omega}$$

Where:
$A$: filter attenuation in dB
$\Delta\omega$: Transition band of filter in radians/sample ($2\pi = $ sampling rate).

The OP's case isn't a simple low pass or high pass filter, but this illuminates that the tighter the features are in frequency, the more taps will be needed.) Similarly the negative derivative of phase with respect to frequency is the group delay of the filter by definition. The OP does not care about phase but for general purposes we see how the slope of the phase vs frequency also drives the number of taps since the filter delay must be within the span of the filter.

The approaches are further detailed below by their common name, but as far as tap estimation, the approach above is reasonably close for the windowing, least squares, and equi-ripple design approaches. For frequency sampling a significantly higher number of taps may be needed  (could be 2x - 3x).

It would be common to start with a much higher number $N$ and then reduced the number of taps once the relative magnitude of the tails at the edge of the filter is reviewed. This can be reduced by pursuing the same algorithm with a smaller $N$, or truncating and windowing the solution achieved with a higher $N$.

Additional details on estimating the number of taps needed:

How many taps does an FIR filter need?

Design Coefficients:

The four common approaches for designing the FIR filter coefficients from target specifications are as follows:

Windowing (OP Option #1)
Frequency Sampling (OP Option #3)
Equiripple (Parks-McClellan, Remez Exchange Algorithm)
Least Squares (OP Option #2)

Windowing: (OP Option #1)

This method is to window the sampled continuous-time impulse response of the filter. The impulse response is determined from the inverse Fourier Transform of the frequency response, which is sampled and delayed to be causal and windowed to select the dominant response features based on the number of taps used.  Just selecting samples of the impulse response with no further modification is multiplying the impulse response by a rectangular window in time. Multiplication in time is convolution in frequency, so the desired frequency response will be convolved with a Sinc function in frequency. The Sinc function has the tightest main lobe for a given time span ($N$ samples), but has relatively high sidelobes. This results in the closest match of transitions but would have more ripple in the passband and stopband compared to specialized windows.

This is demonstrated below with an example of designing a perfect (target design, not achieved) low-pass filter. The impulse response of the low pass filter is a Sinc function, which is delayed in time in order to be causal (resulting in a linear negative phase slope versus frequency) and then truncated with a rectangular window to given the number of taps desired.

Once truncated, the other Sinc function in frequency due to the window (not the Sinc in time due to the desired filter response; in this case slightly confusing since two Sinc functions are involved) will convolve with the desired response, widening the bandwidth slightly and causing ripple in the passband and stopband. This is often improved by simply multiplying the desired impulse response with a different window shape that by design has very low side-lobes. The trade of doing this is a wider main-lobe which will serve to shift the exact locations of frequency transitions but will provide a much better match in areas where the frequency response is not changing rapidly. Increasing the number of taps resolves both cases.

Summary: Windowing is quick and easy if the inverse Fourier Transform of the target response is easily solved, and provides a reasonable result. It is sub-optimal in that other design approaches can achieve closer accuracy to target specifications with less number of taps.

Frequency Sampling: (OP Option #3)

This is the method I mostly discourage anyone from using (although I ultimately used that to demonstrate a quick approach for the OP's more challenging filter at his other post here: How to invert FFT these magnitude equations to get FIR filters?. showing that it does have utility if a quick and useful result is needed).

The benefit of this approach is that it is the simplest: You use the inverse DFT to create the coefficients of the filter. This will result in an exact match at the DFT sample locations, but the resulting response will have significantly higher ripple (deviation) for all frequencies in between.

A hybrid approach that would be very viable in this specific application of the OP's is to zero pad the high frequency portion of the target frequency response (which is the middle of the DFT array) allowing for a much longer duration for the impulse response (eliminating time domain aliasing). The resulting impulse response can then be windowed as in the window approach given above. Optionally the achieved impulse response can be resampled and then truncated/windowed. Basically we would be using a longer Inverse DFT to approximate the Inverse Fourier Transform when the result would be too difficult to compute directly (as in the OP's case).

When deriving the coefficients from the ifft, the frequency response needs to include the "negative frequency" components: for a filter with real coefficients this would be a mirror of the DFT samples from $0$ to $N/2-1$ at $N/2$ to $N-1$. Further, when the magnitude response is of only concern, the following will properly center the impulse response in the span of the filter as depicted in the diagrams above:

coeff = fftshift(ifft(ifftshift(mag_response)));

Here mag_response was generate with the frequency axis going from $-f_s/2$ to $+f_s/2$ so the inner ifftshift translates it to the format expected in the DFT going from $0$ to $f_s$, where $f_s$ is the sampling rate. (corresponding to $0$ to $N-1$). The final fftshift converts the noncausal result in the time domain (the coefficients are centered on $t=0$) to a causal filter by centering the coefficients in the center of the filter span.

Filtering by zeroing out FFT bins is the frequency sampling approach, which is further described here including why it is generally a bad idea (meaning will require a lot more taps than the optimal approaches equiripple and least-squares):

Why is it a bad idea to filter by zeroing out FFT bins?

Equiripple:

This is the go-to design algorithm when the desired accuracy result is constrained by peak limits. The equiripple design approach developed by Thomas Parks and James McClellan in 1972 https://en.wikipedia.org/wiki/Parks%E2%80%93McClellan_filter_design_algorithm uses the Remez-exchange algorithm resulting in an optimum filter design given a minimax constraint (minimize the maximum error). Optimum meaning the solution will be the mimimum number of taps given the peak error design constraint. It can be used to match arbitrary shapes and multi-band filters and is supported by functions in MATLAB (firpm), Octave (remez) and Python (scipy.signal.remez). For a more complicated filter it may have challenges converging which would be my only reason to pursue one of the earlier approaches above.

Least-Squares: (OP Option #2)

This is the go-to design algorithm when the desired accuracy result is constrained by root-mean-square (rms) limits. Given the performance of most applications where I need to use filters (wireless communications) is optimized by rms performance and not peak error, this is my filter design of choice. Like the equiripple design, this is an optimal filter design in that the solution will be the mimimum number of taps given the rms design constraint. Unless I was effectively given a "though shall not cross this line" constraint in the filter specifications, I would use least-squares.  It can be used to match arbitrary shapes and multi-band filters and is supported by functions in MATLAB/Octave (firls), and Python (scipy.signal.firls).  For a more complicated filter it may have challenges converging which would be my only reason to pursue one of the non-optimal design approaches above.

Further posts of interest:

FIR Filter Design: Window vs Parks McClellan and Least Squares


Equalize Speaker using Frequency Response Curve

Oops I misread your question and thought you have a microphone recording of the speaker output with given audio, and answered that. Well here is anyway, as some parts still apply:

Like others have suggested in comments, designing a filter is the way to go. I would recommend a windowed FIR filter design rather than Parks-McClellan. From MATLAB documentation:

  firpm exhibits discontinuities at the head and tail of its impulse
  response due to this equiripple nature.

Our ears analyze sound with kind of wavelets so they may detect those discontinuities before and after impulsive sounds. With windowed filter design however you can choose a window that smoothly starts from and ends at value 0, avoiding creating such discontinuities.

Doing division (inverse of the amplitude) is a bit dangerous, because if you divide something by a number close to zero you get a huge value. Discrete Fourier transform results often contain values close to zero. If you have a long audio file then you can split the audio into shorter pieces with say 50 % overlap, multiply each with a window function, and calculate the mean squared magnitude of their discrete Fourier transforms (discarding phase). The square root of this is the root mean square spectral power, which is unlikely to be zero anywhere. Do this for both the speaker and raw audio and bin-wise divide the root-mean-square spectral power of the latter by the first. Still, you may get very large values if the speaker is unable to reproduce very low or very high frequencies. You can set a (soft) limit to how much correction you do at most, and you can even softly ramp the correction to nothing at frequencies you know the speaker simply cannot reproduce. Assuming zero phase for all bins, you can directly inverse discrete Fourier transform the division result to get an unwindowed zero-phase impulse response. You would then multiply the unwindowed impulse response with a window function which can be much shorter than the transform result, giving you a practical, short filter impulse response.

One more trick is to convert the filter to its minimum-phase counterpart, which will completely eliminate pre-ringing (taking place before the main peak in the impulse response). There are specialized methods for doing this that preserve the magnitude frequency response. I would not worry about post-ringing because real-world acoustics will give you plenty of post-ringing and echos anyhow and we are used to masking them.


How to calculate performance of your object recognition algorithm and compare results with others?

The problem seems to be that the authors of the other papers did not include enough information for you (or anyone) to be able to reproduce their work. This could have been oversight, or the editor (or reviewers) may have removed it prior to publication. I'd suggest contacting the authors of the problematic paper and see if they'll answer your question. Be polite (it's always nice to know someone read your paper, so let them know that you have!)


Sounds like you have read reference 8 and understand that for an infinite $N$ you can one could theoretically achieve the CRLB. But clearly in practice you need to drop that idea and think about what your constraints are for your problem. That is, do you really need an estimate with such low variance and if you do then how long are you willing to wait? Since increasing $N$ means you need more and more samples of your signal to come but if you are doing offline processing then maybe you can afford this and can take a very large $N$.


does aliasing occur always if i sample a vibration in real world applications?

That is why, before sampling, a (steep) lowpass filter with cutoff frequency $f_c \leq \frac{f_s}{2}$ shall be applied. Thus, the amount of aliasing will be insignificant.