It looks like you're getting your frequency in Hertz confused with in radians/sec since you have the factor of $2\pi$ in both your sampling period dt and your signal. I rewrote a bit of your code to clarify what I think you really want.

Amp = 1;
freqHz = 12000;
fsHz = 65536;
dt = 1/fsHz;
t = 0:dt:1-dt;
sine = Amp * sin(2*pi*freqHz*t);

N = 65536;
transform = fft(sine,N)/N;
magTransform = abs(transform);

faxis = linspace(-fsHz/2,fsHz/2,N);
plot(faxis/1000,fftshift(magTransform));
axis([-40 40 0 0.6])
xlabel('Frequency (KHz)')

If your sampling frequency is 65536 samples/second, and you want for example a tone at 12 KHz, you can create it as shown. So here your sample period is 1/65536 seconds.

Your expectation of getting two spikes with an amplitude of 0.5 each was correct -- just your generated tone was not.

As for scaling the x-axis to be in Hertz, just create a vector with the same number of points as your FFT result and with a linear increment from $-fs/2$ to $+fs/2$. Note also the fftshift I used in the plot. That's because the output of Matlab's FFT function goes linearly from 0 to fs. I find it easier to visualize having DC centered but either way is fine. Without the fftshift the faxis vector would go from 0 to fs.


Downsampling and spectral stretching is earlier to show in 1D. I'll also analyze with the Fourier transform and discrete-time Fourier transform, not a DFT. The results are the same except for aliasing in time.

To downsample by an integral multiple is to increase the sample period of a continuous-time signal $x(t)$ by an integer $M : T → M T$ . The replication period in the frequency domain is reduced by the same multiple. Rather than the first replica being centered at $\omega = 2\pi/T$ , the first replica of the down-sampled spectrum is centered at $\omega = 2\pi/MT$. Moreover, the amplitude is scaled by $1/M$. The figure illustrates $M = 2$. Sequencing the original series converts $x(nT ) → x[n]$. To downsample the sequence is to introduce a stride $M$ when reading off the values:

$$
x^{\downarrow} [n] = x[nM ]
$$

The effect on the angular sample frequency can be shown by a units analysis,

$$
\left(\frac{φ^↓\; \text{radians}}{1\; \text{sample}}\right)
=
\left(\frac{ω\; \text{radians}}{1\; \text{sec}}\right)
\left(\frac{T\; \text{sec}}{1\; \text{step}}\right)
\left(\frac{M\; \text{step}}{1\; \text{sample}}\right)
$$

or $φ^{\downarrow} = \omega M T$ . The range $ω ∈ [0, 2π/M T ]$ corresponds to $φ^↓ ∈ [0, 2π]$, thus $φ^↓$ is a stretched version of $ω$. The figure below compares $x[n]$ to $x^↓ [n]$ and $φ$ to $φ^↓$ . $x^↓ [n]$ has fewer samples than $x[n]$, and has a broader spectrum when plotted in $φ$ because the radians per sample has increased by $M$.

So in the 1D analogy to your case, refer to the continuous-frequency spectra $\omega$ above. Downsampling by 2:1 creates a spectral replica centered at $2\pi/MT$, $M=2$. That is what your first lines of code do (in 2D). Then, the spectrum is stretched so that the periodicity of the spectrum is restored to $2\pi$. This is shown to the right.


FDLS is similar to rational polynomial interpolation in that the interpolation quickly becomes unstable if the target polynomial degree is set too high.  (e.g. just a tiny amount of variation or noise in the spec at a sub-optimal grid spacing can cause extremely wild excursions in the graph of a high order polynomial).  Thus it is usually better if the FDLS order you choose is overdetermined by your specification points to some degree (usually selected by trial-and-error).

The jagged-ness of a frequency response is more related to how close a group of poles or zeros are to the unit circle and each other, rather than to their absolute number.


Undesirable zeros in sampled signal using train of impulses

If you just dont want to see them in your figure, just replace them with nans. Lets say that your signal is s. Then

s(s==0)=nan;
plot(s);


Camera calibration/pin hole camera model and working out 3d position

If you have the extrinsics then it is very easy. Having extrinsics is the same as having "camera pose" and the same as having the homography. Check this post in stackoverflow.

You have extrinsics, also called camera pose, which is described as a translation and a rotation:

$\displaystyle Pose =\begin{bmatrix}R|t \end{bmatrix} = \begin{bmatrix}R_{11} &R_{12}&R_{13}&t_x\\R_{21}&R_{22}&R_{23}&t_y\\R_{31}&R_{32}&R_{33}&t_z \end{bmatrix} $

You can get Homography from Pose this way:

$\displaystyle   H = \frac{1}{t_z}\begin{bmatrix}{R_{1x}}&{R_{2x}}&{t_x}\\{R_{1y}}&{R_{2y}}&{t_y}\\{R_{1z}}&{R_{2z}}&{t_z}\end{bmatrix}$

Then you can project your 2D points into the corresponding 3D points by multiplying the Homography by the points:

$p_{2D}=\begin{bmatrix}x &y &1\end{bmatrix}\quad$ add $\quad z=1\quad$ to make them homogeneous

$p_{3D}=H*p_{2D}  $

$p= p / p(z)\quad$   Normalize the points


Which order to perform downsampling and filtering?

You need to filter first and then downsample. Otherwise, you will run into aliasing problems. I.e. frequencies that are above 30 Hz will create images within your frequencies of interest. You can consider the little script below to compare both methods:

Fs = 128.0
t = np.arange(0, 10, 1/Fs)
signal = np.sin(2*np.pi*10*t) + np.sin(2*np.pi*50*t)
sigma2 = 0.5

rx = signal + np.sqrt(sigma2) * np.random.randn(len(signal))

downsampled = rx[::2]
(b1, a1) = scipy.signal.butter(6, 30/(Fs/4))
down_and_filtered = scipy.signal.lfilter(b1, a1, downsampled)

(b2, a2) = scipy.signal.butter(6, 30/(Fs/2))
filtered_and_down = scipy.signal.lfilter(b2, a2, rx)[::2]
t2 = t[::2]

plt.figure(figsize=(20,20))
plt.subplot(221)
plt.plot(t, rx)
plt.plot(t, signal)
plt.xlim((1, 1.3))
plt.title('Received signal')

plt.subplot(222)
plt.plot(t2, down_and_filtered)
plt.plot(t2, filtered_and_down)
plt.title('Compare of both methods in time domain')
plt.xlim((1,1.3))

f = np.linspace(-Fs/2, Fs/2, 4*len(t))
plt.subplot(223)
plt.plot(f, np.fft.fftshift(abs(np.fft.fft(rx, 4*len(t)))))
plt.title('RX spectrum. Note the two peaks, one is interference, having frequency over 30Hz')

f2 = np.linspace(-Fs/4, Fs/4, 4*len(t2))
plt.subplot(224)
plt.plot(f2, np.fft.fftshift(abs(np.fft.fft(down_and_filtered, 4*len(t2)))), '-o')
plt.plot(f2, np.fft.fftshift(abs(np.fft.fft(filtered_and_down, 4*len(t2)))))
plt.title('Spectrum of both alternatives. NOte the blue curve has a wrong component')

In this script you have a signal that is composed of two sines of frequency 10 and 50Hz plus some noise. In the FFT of the full signal, both spectral lines clearly occur. Now, what you want to have after your signal processing is only one spectral line at 10Hz (because the 50Hz signal should be filtered out). HOwever, if you downsample first, the 50Hz wave is mirrored to $14Hz=64Hz-50Hz$ and cant be filtered out subsequently. Hence, you need to do the filtering before and then downsample your signal.


Deconvolution of signal which is broadened by laser spot size

There was a similar question here a while back with a conceptually similar problem.
Link to the earlier question.

That one was about scanning a target with an electron beam.
The question there was how to determine the true shape of the beam if you use a known target.  You would seem to have a beam with a known shape.
You've got it a bit easier, since you are dealing with a single dimension rather than a 2D shape.

The easiest way to start would be to look at examples for the Matlab deconv method.

Basically, you need two vectors - one that represents the beam, and one that represents the object as scanned with the beam.  Feed them into the deconv method, and you get a vector that represents just the object.

One way to go at it would be to create a vector for an object of known size - say you have a sample of exactly 5 microns width.  You make a vector representing say 10 microns, and you put a 5 micron wide step in the middle of it.  Now you scan for real a 5 micron wide object, and deconvolve that vector with your generated 5 micron vector.  Now, you have the beam vector to use when measuring unknown objects.

If you don't have Matlab, you could use Gnu Octave which works almost just like Matlab - many Matlab examples Just Work in Octave.

Example Octave code:

calibration_target=[1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1];
simulated_beam=[1,1,2,2,2,1,1];
calibration_scan=conv(simulated_beam,calibration_target,"full");
[recovered_beam_shape,beam_remainder]=deconv(calibration_scan,calibration_target);

test_item=[0,0,0,1,1,0,0,0];
simulated_testscan=conv(simulated_beam,test_item,"full");
corrected_test_item_scan=deconv(simulated_testscan,recovered_beam_shape);

subplot(4,2,1);
plot(calibration_target);
title("calibration target");

subplot(4,2,2);
plot(calibration_scan);
title("simulated scan of target with simulated beam");

subplot(4,2,3);
plot(simulated_beam);
title("simulated beam");

subplot(4,2,4);
plot(recovered_beam_shape);
title("recovered beam shape");

subplot(4,2,5);
plot(simulated_testscan);
title("simulated testscan");

subplot(4,2,6);
plot(corrected_test_item_scan);
title("corrected test item scan");

subplot(4,2,7);
plot(test_item);
title(" test item true shape");

Result:

Added to answer questions from the comments.

Convolution and deconvolution work on uniformly spaced vectors.  All you need are the intensity values - drop your x coordinates.  If your data points aren't uniformly spaced, then you've got a problem.  If you interpolate your points to get uniform spacing, then you will probably introduce some fuzziness.  I don't know how much of a bad effect that would have on the deconvolution.  As for determining the vector to correct for the beam, I would scan a target with a known shape and sharp edges.  Then I would deconvolve the scanned data with an ideal image of the target - basically a vector as wide as the scanned object, with values of zero for the areas outside of the target, and values of 1 inside the target.  For example, a target 5 microns wide scanned with 1 micron stept would be [0,0,0,0,0,1,1,1,1,1,0,0,0,0,0]  Since you did a scan with a known object in order to generate your gaussian, I assume this is possible.
The "size" shouldn't matter - assuming you mean the amplitude or intensity.  I would scale the generated data to make 1 represent the scanned area at maximum intensity.  If you meant the number of elements, then that would matter.  You need to overscan by enough to be sure that your beam actually leaves the target to each side.  If you don't overscan enough, then you won't be able to get a true image of your target.


If you want to vertically upsample an image, one method possibility better than linear interpolation is to FFT, zero pad, and IFFT; or a similar process using the DCT/IDCT and/or sub-blocks of the image.


Converting ADC output to fixed-point format

You don't have anything to do, because $\rm Qm.n$ is just an "interpretation" of an integer. The only thing that really matters is that you need to correctly apply a shift operation when multiplying numbers; or when adding/subtracting numbers with different values of $\rm (m, n)$ - You can think of the $\rm Qm.n$ representation as a "tag" you attach to number, a bit like physical units.

In your case, the output of your ADC is an integer between $-32768$ and $32767$. This number can be interpreted as a $\rm Q15$ fixed-point number between $-1.0$ and $0.99997$. Or as a $\rm Q2.14$ fixed-point number between $-4.0$ and $3.99987$. Or as a $\rm Q8.8$ number between $-256.0$ and $255.992$. You don't have to do any arithmetic operation.


See related question: Is initializing a digital filter's output with no "momentum" a non-trivial task?

The easy way is to subtract, prior to filtering, from each input sample the first input sample. This will eliminate the step in the beginning, avoiding the wavy step response. That's roughly what the other original poster said they did. Or because you have noise, you can take the average of the first few input samples as a better estimate of the DC bias and subtract that from each input sample, prior to filtering.

If you want to save the work of subtracting then changing the filter's state variables in the beginning is the way to go, and my answer to the linked question shows how to do that (with $x$ as the estimated DC bias and $y=0$ as you have a high-pass filter with zero DC response but this can be validated). Not sure if you can easily touch the state variables in your prototyping environment without rewriting the filter code.


Determining Frequency and Period of a wave

Interesting project you have going on there! :-)

From a signal analysis POV, this is actually a simple question - and yes, you are right that you would utilize the FFT for this frequency estimation problem.

I am not familiar with R, but what you essentially want to do is take the FFT of your temperature signal. Since your signal is real, you will get a complex result as the output of your FFT. You can now take the absolute magnitude squared of your complex FFT result, (since we dont care about phase). That is, $real^2 + imag^2$. This is your power spectral density.

Then, very simply, find the max of where your PSD sits. The abscissa of that max will correspond to your frequency.

Caveat Emptor, I am giving you a general outlook, and I suspect the result of the FFT in R will be normalized frequency, in which case you would have to know your sample rate, (which you do), in order to convert it back into Hz. There are many other important details that I am leaving out, such as your frequency resolution, FFT size, and the fact that you probably want to de-mean your signal first, but it will be good to see a plot first.

EDIT:

Let us take your signal into account. After I de-mean it, it looks like this:

You want to de-mean because the DC bias is a frequency of 0 Hz technically, and you do not want it to drown out other frequencies of interest. Let us call this de-meaned signal $x[n]$.

Now, when you take the FFT, you must take heed of some details. I use an FFT length of 10 times the length of the signal, so we can say that $N_{fft} = 10 * 3600 = 36000.$ Any FFT size set greater than the size of the original signal has the effect of interpolating your FFT result in the frequency domain. While this doesnt add any information for you from an information theory POV, it does help you better discern where exactly your true peak lies, especially in cases where it lies between two frequency bins. To get true better frequency resolution, you want to get more data. (i.e, let the sensor run for a longer period of time).

Moving along, from the text file, I see that your sensor is taking temperature readings every 2 seconds. This means that your sampling rate, or $f_s = 0.5 Hz$

Thus, let us take an FFT of $x[n]$, and call that result $X(f)$. This result will be complex, and in fact, it is conjugate symmetric, but thats not important here. (It just means for your purposes, half of it is redundant). Now if we take $10log_{10}(|X(f)|^{2})$, this is the power spectral density (in log scale). The result looks like so:

You can see how it is symmetric. If you ignore the last half, and just look at the first half and zoom in you, can see this:

So you have a peak at index 70! But what is index 70 in real frequency terms? Here is where you want your frequency resolution. To compute that, we simply take $\frac{F_s}{N_{fft}} = 1.3889e-005 Hz$. This simply means that every index represents a frequency that is a multiple of this number. Index 0 is $0 * 1.3889e-005 = 0 Hz$. Index 1 is $1 * 1.3889e-005 = 1.3889e-005 Hz$. Index 70 is $70 * 1.3889e-005 = 9.7222e-004 Hz$.

We are almost done. Now that you have this figure, we can convert it to something more palatable. This figure is just saying that your signal completes 9.7222e-004 cycles every second. So we can ask, how many minutes does it take to compute one cycle? Thus, $\frac{1}{(9.7222e-004)*60} = 17.14$ minutes per cycle.


Need information about following channel model

The channel impulse response represent the weighted gains and phases at different delays, so yes this is a multipath channel. You can simulate the channel (in Matlab as you requested, what I show below was done in Octave with compatible results) by first creating the complex channel weights:

channel = [c1, c2, c3, c4, c5, c6, c7];

where $c1 = -.005 - j0.004$, $c2 = 0.009 + j0.030$, etc

The following are examples of inspecting and using the channel:

filter

Use the filter command to determine the complex baseband received signal from the complex baseband transmitted signal (analytic signal):

rx = filter(channel, 7, tx)

Where rx is the received signal and tx is the transmitted signal. The denominator factor of 7 is somewhat arbitrary and is just a scaling.

freqz

Use the freqz command to plot the magnitude and phase frequency response of the channel.

freqz(channel)

The horizontal axis is normalized angular frequency from DC to $2\pi$ where $2\pi$ represents the sampling rate used. (So if the tx data was sampled at 1 MHz for example, the right side of the plot would be 1 MHz). Because it is a complex spectrum the full frequency range from DC to the sampling rate is unique as shown.

grpdelay

Use the grpdelay command to determine the group delay of the channel, which is the delay of the signal in time or in samples at each frequency. In the phase vs frequency plot above it appears the channel may be linear phase, however a group delay plot (group delay is $d\phi/d\omega$) would give us a better indication since it will remove the linear trend.

grpdelay(channel, 1, 512, 'whole');
xlabel('Normalized Frequency (rad/sample)')
ylabel('Delay (samples)')
title('Group Delay')

From the above results we see that the channel distortion is not significant or challenging both in magnitude (+/- a few dB) and group delay (+/- half a sample), in comparison to a highly frequency selective channel that would have deep nulls throughout the passband magnitude response. The delay distortion is more pronounced if a lower sampling rate is used, so depends on how many samples per symbol there are (The delay can be converted to seconds by dividing by the sampling rate).

See this post on implementing an LMS Equalizer to compensate for this channel in the receiver (using the Wiener-Hopf equations): Compensating Loudspeaker frequency response in an audio signal .

For "blind" equalization using this approach, I would use receiver decisions to recreate the expected transmitted signal (including any pulse shaping etc). I did not read your reference link but you could do that as a comparison if that wasn't the technique used directly. A comparison should include performance under lower SNRs.


Is zero forcing Optimization technique or Equalization technique?

Zero-forcing is not an optimiziation technique and neither an equalization technique. It is simply one of several possible criteria to design (or optimize) an equalizer. A zero-forcing equalizer is designed to completely remove intersymbol interference, but does not take noise into account. Another, more practical criterion is the MSE (mean-square error) criterion, which results in an equalizer that minimizes intersymbol interference and reduces noise power at the same time.

I would suggest you consult one of the many good books on digital communication to learn more about this quite complex topic. Here are just two which I like:

Digital Communication by Lee and Messerschmitt
Digital Communications by Proakis


Alternative to BIBO stability of a system

For linear systems, BIBO stability is the most useful and practical criterion. For systems described by rational transfer functions it coincides with the condition that all transfer function poles must be located in the left half of the $s$-plane (for continuous-time systems), or inside the unit circle of the $z$-plane (for discrete-time systems). For linear systems, I do not know of any other stability criterion that makes any sense and does not coincide with the BIBO-stability criterion. Since most basic DSP texts focus on linear (and especially time-invariant) systems, you will only find the notion of BIBO-stability there.

Things look a bit different for non-linear systems. There you have several reasonable stability criteria that do not all lead to the same basic criterion as is the case for linear systems. One important notion for non-linear systems is input-to-state stability, which basically means that for zero input, the system is stable about its zero state, and that well-behaved and bounded input signals produce a bounded state trajectory. This article reviews some of these concepts.

But if you're mainly interested in linear systems, BIBO-stability is all you need.


How to take the FFT of an undersampled signal

If there is other frequency content in other Nyquist regions, you will see those alias into your digitized signal also; there's no free lunch here. If you wish to use bandpass sampling, then you need to first apply a bandpass anti-aliasing filter to suppress the unwanted content.


Duality in DFT would mean that if $x[n]$ has DFT coefficients as $X[k]$, then DFT of $X[n]$ would be $Nx[(N-k) \mod N]$

Proof:
Given, $$X[k] = \sum^{N-1}_{n=0}x[n]e^{-j\frac{2\pi}{N}nk}, k=0,1,2,3,...,(N-1)$$
If we take DFT of the sequence $X[n]$, then what we get is the following :
$$Y[k] = \sum^{N-1}_{n=0}X[n]e^{-j\frac{2\pi}{N}nk} = N \left(\frac{1}{N}\sum^{N-1}_{n=0}X[n]e^{j\frac{2\pi}{N}n(-k)} \right)$$
Notice that the expression between "()" is the synthesis expression with $k^{th}$ frequency index replaced by $(N-k \mod N)^{th}$ index. Because, in DFT expression, $k$ can only take indices $0,1,2,3,...,(N-1)$, hence, we cannot have $-k$ as frequency index, but rather, $((N-k) \mod N)$
$$Y[k] = N \left(\frac{1}{N}\sum^{N-1}_{n=0}X[n]e^{j\frac{2\pi}{N}n(N-k)} \right) = Nx[(N-k) \mod N]$$
So, the way to interpret this is, you get a scaled and inverted sequence back when you take DFT of DFT, but $x[0]$ remains at $0^{th}$ index.
So, DFT of DFT of time-domain $x[n]$ gives $\{Nx[0], Nx[N-1], Nx[N-2], ..., Nx[2], Nx[1]\}$
Scaling by $N$ is the consequence of not dividing by $\frac{1}{\sqrt{N}}$ when taking DFT and incorporating this factor into IDFT expression.


Conflict with the properties of convolution

I believe that the problem you have with this exercise is mainly caused by bad notation. Writing convolution as

$$z[n]=x[n]\star y[n]\tag{1}$$

(which I think is bad notation) it is tempting to believe that

$$z[n-1]\stackrel{?}{=}x[n-1]\star y[n-1]\tag{2}$$

which is not the case.

If you take the definition of (discrete) convolution

$$z[n]=\sum_kx[k]y[n-k]\tag{3}$$

then it's obvious that

$$z[n-1]=\sum_kx[k]y[n-1-k]=\sum_kx[n-1-k]y[k]\tag{4}$$

i.e., only one of the two signals needs to be delayed by one sample, not both.

A better notation is

$$z[n]=(x\star y)[n]\tag{5}$$

Now you can define $x_1[n]=x[n-1]$ and $y_1[n]=y[n-1]$ and you have

$$z[n-1]=(x\star y_1)[n]=(x_1\star y)[n]\tag{5}$$

If you view $z[n]$ as the output of a system with impulse response $y[n]$ and input $x[n]$, then it should also be obvious that delaying $z[n]$ by one sample is equivalently achieved either by delaying the input signal, or by adding a delay to the impulse response, but not by doing both.


Why idft(dft(a) * dft(b)) not equal to convolve(a, b)?

ifft(fft(a) * fft(b)) performs a cyclic convolution, convolve apparently zero-pads the inputs. If you pad both arrays with zeros, the result should be the same:

a = [0,0,0,1+0j, 2+0j,0,0,0]
b = [0,0,0,4+0j, 5+0j,0,0,0]

print list(ifft(fft(a) * fft(b)))
print list(convolve(a, b))


whether the system is linear or not for the given problem

How can we calculate output at 0 input if the system depends on future or past or both?

Well, surely zero input just means;
$$
x(t) = 0~~~~\forall t
$$
and the $\forall t$ means for all time: positive and negative.

Substituting that into the equation:
$$
y(t) = x(t+1) + x(t-1) = 0 + 0 = 0
$$

So the system is homogeneous.

Well, as @Dilip points out, this isn't sufficient for homogeneity: we need the output to be $a y(t)$ for all inputs $a x(t)$. The case above just looks at $a=0$.

The next question is: does it satisfy additivity?

If
$$
x_{\tt total}(t) = x_1(t) + x_2(t)
$$
then
$$
y_{\tt total}(t)  = x_{\tt total}(t+1) + x_{\tt total}(t-1)\\
= x_1(t+1) + x_2(t+1) + x_1(t-1) + x_2(t-1)\\
= y_1(t) + y_2(t)
$$
where
$$
y_1(t) = x_1(t+1) + x_1(t-1)\\
y_2(t) = x_2(t+1) + x_2(t-1)
$$
so it satisfies additivity also.


Given that you know the autocorrelation function, you can use the Yule-Walker equations for a generic AR model of order $p$. As you are interested in the case $p=2$, there's no need to complicate things that much.

As stated in this paper, the equations for the case of an AR(2) are given by:

$$a=\frac{r_1(1-r_2)}{1-r_1^2}$$

$$b=\frac{r_2-r_1^2}{1-r_1^2}$$

where $r_i$ are the coefficients corresponding to the autocorrelation function.


reference : Digital Image processing - Rafael Gonzalez

Steps :
The whole of filtering process can be summarized in the following points :
1> Given an input image f(x,y) of size M X N, obtain the padding parameters P and Q as P =2M and Q = 2N

The problem with using DFT is that both the input and output images are periodic. This causes interference between adjacent periods if the periods are close with respect to the duration of the non-zero parts of the functions. So we padd  the image by using a padding function P having zeros equal to twice the image length.

2> Form a padded image f’(x,y) of size P X Q by appending the necessary numbers of zeros to f(x,y).

3> Multiply f’(x,y) by (-1)^x+y to center its transform

4> Compute the DFT, F(u,v) of the image.

f = imread('pic.jpg');
figure(1);
imshow(f);
title('original image');

PQ = paddedsize(size(f));
[U,V] = dftuv(PQ(1),PQ(2));
 F = fftshift(fft2(f,PQ(1),PQ(2)));

 S =  abs(F);
 ag = angle(F);

[l,p] = size(S);
d = max(max(S));
c = min(min(S));

for l = 1:1:l
    for p = 1:1:p
        S1(l,p) = (S(l,p)-c)*(255/(d-c));
    end
end

figure(2);
subplot(2,2,1);
imshow(uint8(S1));
title('Magnitude response without log transform');

% since dc dominates the values of spectrume, the dynamic range of other intensities in displayed image is compression. to bring out those details we perform log transformation.

[m,n] = size(ag);
d = max(max(ag));
c = min(min(ag));

for l = 1:1:l
    for p = 1:1:p
        ag1(l,p) = (ag(l,p)-c)*(255/(d-c));
    end
end

subplot(2,2,2);
imshow(uint8(ag1));
title('Phase response without log transform');


Your low pass should be higher - probably 40 or so.
Your high pass can be higher - around thirty
Your gain should be variable else you will amplify it into distortion - what you want is called normalisation.  You will need to find the maximum volume of your signal and calculate a scaling factor that will bring this up to full scale.  (Say full scale is 1.0 and your maximum value is 0.01 then multiply all samples by 100.)

When I experiment with heart sounds, I find the beat itself to be around 34Hz - I usually do a band pass from 30 to 40 in order to isolate just the beat.

If you do this, you will have a clean recording of ONLY the heartbeat.  If all you want is to count heartbeats or listen to them you will be fine.  If you need to listen for heart defects or other things, then you will need to change the parameters.  You will need higher frequencies and a wider bandwidth for those types of things.

One other thing you need to consider:  30 to 40 Hz doesn't really come through headphones and ears all that well.  Even when you isolate the hearbeat and amplify it to full scale, it is still hard to hear.  Earphones will attenuate it some, and your ears even more so.

As to actually doing the filtering:

You can use IIR filters at higher orders to isolate the heart sounds.  Avoid using chebyshev filters, though.  The "ring" badly and can make one beat "smear" into the next.
If you use an FIR, you may want to downsample first.  It depends on your sampling rate.  If your sampling rate is high (say the typical 44100) then an FIR than can bandpass something around 30Hz will have to be VERY long.  If you downsample to 1000 Hz, then shorter filters can be used.  You could get away with maybe 100 taps at 1000Hz.  Downsampling consists of doing a low pass filter to remove everything above half of your new sampling rate, then "throwing away" a lot of your samples - downsampling by two would throw away every second sample, downsampling by 4 keeps one out of four samples.

When I do audio in Java, I usually import libpd, and use pd to develop the actually processing stuff.  Then the same pd "patch" can be used with libpd in the java programm.
pd  libpd

The order in which you apply the filters (for your  bandpass) shouldn't matter.


What is the time it takes to measure the SINR in LTE?

So you're giving your answer yourself in the question:

  In general, how the standard LTE works to evaluate the SINR at the receiver?

The point is that measures like SNR and SINR are dependent on what you consider your signal of interest, what you consider interference and what you consider noise¹.

LTE actually specifies multiple different things that describe "received signal strength", whcih you can pu in relation to unwanted effects, and all of these could be labeled as SINR (but have different values for the same transmission!). The two "S" in ("SINR") measurements you probably care about are

RSRP, Reference Signal Received Power: Simply put, LTE uses distributed pilot symbols to figure out the channel state; just adding up the strengths of these will give you the magnitude of the impulse response over all OFDM subchannels.
RSRQ, Reference Signal Received Quality: The reference tone powers above don't say overly much about how well-suited the channel actually is – they are really just a sum. But in an OFDM system, where all but let's say 5% of subchannels are heavily attenuated, but these 5% are extremely well power-conserving, you just won't get the same performance in a system where all channels are attenuated, but things are more evenly distributed, even if the two systems have the same RSRP; hence, RSRQ takes other things into consideration, too.

You can (weightedly) subtract these from RSSI, which is also specificially defined for LTE, to get an estimate for the "IN" part of "SINR".

Frankly, RSRQ and RSRP are probably much better suited to describe a complex multicarrier system than the simplistic SINR – so if I could, I'd just try to argue based on these, and drop the SINR terminology.

¹ This is actually the answer I most often give at my job, "SNR depends on what you consider to be signal and noise. No general formula can be given.".


So, you have a bunch of datapoints of the form (x,y), and considering all of those datapoints together, you have the vectors $x$ and $y$.

To do a curve fit, you would like to solve the equation $Aw=y$ for the column vector $w$, which holds the coefficients of your curve fit polynomial.  These coefficients are also the weights of the basis vectors in the matrix $A$.

The matrix $A$, for an $n^{th}$ order approximation, will hold $n+1$ columns and have a number of rows equal to the length of $x$ and/or $y$ (the number of datapoints).

The first column of $A$ is the vector $x$ raised element-wise to the power of 0. So, it's just a vector of 1's.  This is the basis associated with the offset or y-intercept.  The second column of $A$ is the vector $x$ raised element-wise to the power of 1. So, it's just the vector $x$.  The third column of $A$ is the vector $x$ raised element-wise to the power of 2.  You see the pattern?

To solve for the weight vector, $w=pinv(A)y$, the pseudoinverse of $A$ times $y$. Note, this doesn't actually solve the equation $Aw=y$, but performs an optimal projection of your data onto the vector space described by $A$.

EDIT:  $pinv(A)=(A^TA)^{-1}A^T$


Books/resources for implementing various mathematical functions in fixed point arithmetic for DSP purposes


Differences Using Maximum Likelihood or Maximum a Posteriori for Deconvolution / Deblur

You can think of MAP as a regularization of the ML.
Just like you have regularization for Least Squares Problem (They can be built, mostly, as MAP problem).

The nice thing is that, as always, the best regularization is more data, namely, in most case when there is a lot of data they collide (Namely, low sensitivity fir the Posterior PDF).

So they differ mainly when there is (Relatively) low amount of data.
Now, when you have low amount of data and you know nothing about the parameters you're trying to estimate, ML is the way to go.
If you have low amount of data yet some prior knowledge about it, or reasonable assumption, make those assumption as regularization.
It even good to make the effort and describe this knowledge in Posterior PDF form.

Update

Found really nice tutorial about the topic - The Truth About Priors and Overfitting.


Coherent averaging of polyphase components of a signal

Yes. Simply sum all the polyphase outputs and the sum result will have higher resolution.

Consider that each polyphase output is a delayed version of the same signal, so that if you commutated through all the outputs, you would get a higher sampled version of your same signal and the quantization noise of this signal would be approximately white across this higher digital frequency band. If you passed this signal through a low pass filter prior to taking every nth sample in a decimation process you would indeed benefit by an increased dynamic range to the extent the quantization noise components in each sample are independent. (they should be; the delayed replicas at each output would have independent quantization noise and coherent signal components). So the decimation filter could be a simple moving average in which case it is easy to see that your best solution is to sum average the outputs of your polyphase filters and use that average as the final result rather than just selecting any one of them.  The moving average introduces passband droop- so with more complication you could do more elaborate filtering. I would be inclined to use a CIC filtering approach as the moving average followed by a simple 3 tap inverse sinc shaper if I was concerned about the resulting passband droop.


High-pass filter with quantized (fixed-point) coefficients

One problem with your design is that the quantized numerator coefficients don't add up to zero anymore, so you lose the desired notch at DC. That can be done differently, because even with quantized coefficients you can have a (double) zero at DC. Just make sure that $b[0]=b[2]$ and $b[1]=-2b[0]$.
The real problem is the denominator polynomial. It is well known that the available pole locations resulting from quantizing the polynomial coefficients are very sparse at low frequencies.
Note that the maximum magnitude of the ideal coefficients is close to $2$, so in order to quantize the coefficients to $16$ bits, you need to scale by $2^{14}$ instead of by $2^{15}$. This will result in the correct value range of the quantized integer coefficients:

bb = round( b * 2^14 ) * 2^(-14);
aa = round( a * 2^14 ) * 2^(-14);

The resulting quantized numerator coefficients guarantee a double zero at DC. However, the quantized denominator coefficients result in two real-valued poles, one of which is on the unit circle. I.e., there is a pole-zero cancellation and you end up with a first-order filter.
So one of the options you have is to use a first-order filter instead of a second-order filter. The numerator coefficients are trivial in this case:

b1 = [1,-1];

The denominator coefficient can be chosen experimentally using some negative value close to (but smaller in magnitude than) $-1$. Of course that value must be representable using $16$ bits:

a1 = [1, round( -r * 2^(15) ) * 2^(-15)];

The plot below shows the frequency response of the ideal second-order filter (blue) and of the quantized first-order filter with $r=0.998$. You have to judge for yourself if the deviation from the ideal response is tolerable for you application or not.

Another solution is to use a different filter structure which is more robust to coefficient quantization (at least at low frequencies). Such a structure is the coupled-form structure, which is described in most textbooks on DSP (e.g., in Oppenheim's Discrete-Time Signal Processing).


Why can Fourier series be used on a non-repetitive function?

Fourier series can be calculated for any signal in close section (from $T_0$ to $T_1$). In fact, I was taught that this is how Fourier series is defined.

But Fourier series is made up of ($T_0 - T_1$) cyclic functions, therefore, if you go beyond the original section you will just get repetitions as you displayed.


How do I find transfer functions from a state space representation?

or the transfer function from the input  to the output, how can I do this ?

Is that not what your $G$ is?  If you want to find the transfer function from just one element of $u$ to the output, then either delete the columns of $B$ that don't pertain to that element of $u$ and get your transfer function, or just look at the column of the transfer function that matches the element you want.

  but now, suppose I want to obtain the transfer function from the disturbance to the output,

Then you would make a column for $B$ (or make a new $B$) that represents the effect that the disturbance has on the system, and extract the transfer function from that.

Edit

I neglected to include the actual math that Matlab is doing under the hood.  This is twice bad -- once because I did it, and twice because I really don't like people who just push the Matlab buttons without understanding what's actually going on.

If you have a system in state space representation:
$$\begin{split}x_k = A x_{k-1} + B u_k \\
y_k = C x_{k-1} + D u_k\end{split}$$ then you can take the $z$ transform:

$$\begin{split}X = A\ X(z) z^{-1} + B\ U(z) \\
Y = C\ X(z) z^{-1} + D\ U\end{split}$$
Then (leaving out steps, the grader will have words with me, but it's all linear algebra) you can solve for $X/U$:
$$\frac{X(z)}{U(z)} = C \left(I z - A\right)^{-1} B + D$$

Note that a nice thing about doing it this way is that if $B$ and $C$ happen to be single column and single row, then you get a nice, traditional scalar transfer function -- but if they're multi-dimensional, you get a very natural matrix representation for a transfer function that just drops naturally out of the math.


How can I create a jamming signal for desired frequency?

The aim is to get SINR of the desired receiver low. Take a high power signal and spread it in the desired bandwidth. For ex: take any bit sequence and spread it in a bandwidth using CDMA  like pseudo noise sequence. The bandwidth in which the sequence is spread will depend upon the chipping rate of the pseudo noise sequences. If the power of the interferer is high compared to the signal power (3db or more) then the intended signal will be degraded. Even if the GPS is using some pseudo noise sequence this will still be enough to cause interference as long as interferer power is high. Very similar to the Near Far Problem in CDMA

PS: Hope its only for study purposes :)


1) It depends on the application, but in general, both. In some cases we can whiten the noise, say in the following situation:

$$Y = x+ \eta\text{ where }\eta \sim N(0, \Sigma)$$

where we want to estimate $x \in \mathcal{A}$ for some set $\mathcal{A}$ from $Y$. In this case, we multiply both sides by the inverse of the square root of $\Sigma$, and we have an estimation problem where $\Sigma^{-1/2} x \in \Sigma^{-1/2} \mathcal{A}$.

2) Two random variables are said to be correlated if their covariance is non-zero. In the cause of jointly Gaussian random variables, uncorrelated is equivalent to independent.

3) Uncorrelated things are often easier to deal with since things decouple (for example, when estimating in the MMSE sense, uncorrelated is a form of orthogonality (perpendicularity, geometrically)). Its the next best thing to independence (and in the most important case, the jointly Gaussian case, they coincide). As for when things are Gaussian, one may justify things via central limit theorem or just by saying the answer is usually easy to find, and several estimators coincide in the jointly Gaussian case (MAP, MMSE, MMAE) [In the jointly Gaussian MMSE case, the estimator is the linear MMSE estimator, which is quite easy to calculate and only relies on first and second order moments, which makes life a lot easier instead of relying on higher order moments].

MLE is just a paradigm for non-random parameter estimation (where as in Bayesian settings such as MAP, MMSE, MMAE, the quantity you are trying to estimate is modeled as random). None of these methods impose a constraint of signal+additive noise for parameter estimation.

You should try reading a book or two on the topic.


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

should just be

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

same as the question you asked yesterday


How is it possible to Filter out a person's voice out of 100 of other voices?

If the signal is recorded using just one microphone, you can use methods such as spectral subtraction. This method is more suitable for "constant" noise, like the noise from a fan or an idle engine. Other methods rely on statistics and perceptual models of speech. If the signal is recorded with several microphones, you can use blind source separation for separating the (speech) signals. As it stands today, you won't get perfect results. The typical end-result is always a trade off between "noise" and clarity of the speech signal of interest. More "noise" suppression --> more degradation of the signal of interest.


When can the impulse response become zero?

The notion of sparsity entails that an object, living for instance in an $n$-dimensional space, can be described (in the suitable basis/frame) by a number $k$ of meaningful components (each above a threshold, or whose combination is close enough to the signal)  "much smaller" than $n$.

When talking about filter identification (adaptive or not), the filter is, at the beginning, unknown. To avoid looking inside an infinity of possibilities, one first restricts  the space of solutions to the subspace of length-$n$ or order-$n$ filters. Either for true physical reasons, out of wishful thinking or heuristic motivations, or because  the optimisation problem is not directly tractable (non-smooth, non-convex), one may call for Ockham's razor, or the law of parsimony, one may want to restrict the solution (not the real object) to be $k$-sparse.

Your example could  be sparse per se, if initially the space of solutions could have up to 100 taps. Additionally, it could perhaps be approximated (in practice) by a sparser filter:
$$y[n] = x[n] + 0.6*x[n-3]$$
discarding the two smallest coefficients.

With this introduction, let us go back to your questions:

I would disagree with "FIR filter = Moving Average and IIR = Autoregressive model": an IIR can be AR/MA, some FIR can be rewritten in an AR fashion. The rest of the question is unclear to me: if the model is linear, the impulse response would be some sparse approximation of the actual coefficients.
The actual coefficients are not becoming zero. But, either they form a  filter short enough to be qualified as sparse, or their estimate has very few non-zero taps
It may happen in some echo cancellation problems: either in a room or on delay lines, one expect that the source bounces after some travelling back and forth. If at each bouncing the reflection is scarce, and there is some attenuation, the  "global" echo filter could look like  (two bounces)
$$y[n] = x[n] - 0.25*x[n-99]+ 0.5*x[n-100]+ 0.25*x[n-101]  + 0.125*x[n-199]+ 0.25*x[n-200]+ 0.125*x[n-201]$$
hence a sparse filter, even with order 201. Outside speech processing, we used such priors in seismic processing of remove multiples, approximate models of bouncing waves, assuming that the filters were sparse at boundaries between rock layers, see A primal-dual proximal algorithm for sparse template-based adaptive filtering: Application to seismic multiple removal, 2014, IEEE Transactions on Signal Processing:

  Unveiling meaningful geophysical information from seismic data
  requires to deal with both random and structured “noises”. As their
  amplitude may be greater than signals of interest (primaries),
  additional prior information is especially important in performing
  efficient signal separation. We address here the problem of multiple
  reflections, caused by wave-field bouncing between layers. Since only
  approximate models of these phenomena are available, we propose a
  flexible framework for time-varying adaptive filtering of seismic
  signals, using sparse representations, based on inaccurate templates.
  We recast the joint estimation of adaptive filters and primaries in a
  new convex variational formulation. This approach allows us to
  incorporate plausible knowledge about noise statistics, data sparsity
  and slow filter variation in parsimony-promoting wavelet frames. The
  designed primal-dual algorithm solves a constrained minimization
  problem that alleviates standard regularization issues in finding
  hyperparameters. The approach demonstrates significantly good
  performance in low signal-to-noise ratio conditions, both for
  simulated and real field seismic data.

Note that parsimony, usually resorting to an $\ell_0$ count measure, is not easy to optimize, as non-differentiable, non-convex, etc.


Your goal is to create a signal with 3 times more samples. The first step is to insert 2 zeros between each sample.
Technically, you now have a signal with 3 time more samples, as you require. However most of the samples are zeros. What you need to do next is to interpolate the 2 samples between each 2 nonzero samples. And this is exactly what the low-pass filter does (or, at least, this one way to look at it).
Digital filters introduce delay to the signal by half the filter's length. I believe this is the reason for the time shifting


How can an audio segment report an integrated loudness (LUFS) greater than 0

I did download this annoying video and analyzed the part with a vuvuzela. You should take a note that:

This part is clipped
The audio is stereo

And that pretty much solves the mystery.

If you calculate the loudness for each of the channels separately you will get approximately -1.4 LUFS. Since this is a stereo track you add those levels together. For incoherent sources, you would roughly expect a 3dB of boost. Lo and behold, the total loudness calculated from the stereo track is approximately 1.6 LUFS (+3.2).
If you think about it, this makes perfect sense. What happens if you take one of those clipped channels and play it back through a speaker? You gonna get some SPL. Now, add a second speaker and play the other channel, you would expect to be even louder (3dB louder to be precise). Therefore LUFS accounts for that.


The correct code is

metrics.Beta = sum(sum((I1 - I_1).*(I2 - I_2)))./(sqrt(sum(sum((I1 - I_1).^2))).*(sum(sum((I2 - I_2).^2))));

The reference you have shown is incorrect, reffr paper:

  2013 August Gregorio - ‘A suitable threshold for speckle reduction in
  ultrasound images’,  IEEE transactions on instrumentations and
  measurement


Spectral flatness or Wiener Entropy for AR(1) and AR(2)

I found the answer:

The AR(1) Gaussian process is defined as:

$$ x(n) = a x(n-1) + \sqrt{1-a^2}\epsilon (n) $$

Where $\epsilon (n)$ is zero-mean stationary white Gaussian noise.

It has PSD:
$$\gamma_{xx}(f)=\frac{\sigma^2_{x}(1-a^2)}{1-2a\cos(2\pi f)+a^2} \ \ \  |f|<1/2$$

Where $$\frac{\sigma^2_{x}}{\sigma^2_{\epsilon}}=\frac{1}{1-a^2}$$

The AR(2) is defined as:

$$ x(n) = a_2 x(n-2)+a_1 x(n-1) + \sqrt{G}\epsilon (n) $$

The definition of $\epsilon (n)$ is the same as above.

Its PSD reads:
$$\gamma_{xx}(f)=\frac{\sigma^2_{x}G}{|1-a_1 e^{-j2\pi f}-a_2 e^{-j4\pi f}|^2} \ \ \  |f|<1/2$$

Here:
$$\frac{\sigma^2_{x}}{\sigma^2_{\epsilon}}=\frac{1}{G}$$

According to this this Paper for an AR(p) Gaussian, the Spectral Flatness Measure (SFM) is very easy to compute:

$$ \hbox{SFM=}e^{-2\rho}\  \hbox{ where } \rho=\frac{1}{2}\log_2\left(\frac{\sigma^2_{x}}{\sigma^2_{\epsilon}}\right)$$

Their SFM is respectively:

$$ e^{\log_2(1-a^2)} $$
and
$$e^{\log_2(G)} $$


Whats the optimal window function to use for analyzing real-time data samples?

A full-width rectangular window (a.k.a no window) would best emphasize spectral content near the edges of an FFT aperture, at the cost of artifacts caused by convolution of the spectum with a Sinc function, also called "spectral leakage".  If you can ignore or compensate for these spectral artifacts then there's little reason to window.

But you might be asking the wrong question.  If you want the fastest possible response to some specific frequency band, you should be using a minimum phase filter (or a bank of them, as needed) or a close approximation to such, not a faster block DFT.  Then you can update with each new sample instead of waiting to process a block, and tune the time-frequency response of the filter(s) for each specific band of interest.


There are plenty of devices that can do this. Simply Google portable recorders. They generally all have a microphone built in. Some have line level inputs, and some can also take XLR inputs and provide phantom power for condenser microphones.

If you don't need the XLR inputs, I can recommend the Zoom H1. It has a stereo microphone, can also take a line level input and is fairly cheap. I have used this to record various foley sounds including bird songs. The sounds save onto an SD card in high quality 24-bit/96 kHz WAV files for analysis on a computer later.

It's best to get a little tripod for it too to stop any handling noise.


Understanding the FFT phase spectrum with a simple example

For sinusoids that are not exactly integer periodic in the FFT length, an FFT measures the phase at a circular discontinuity.  And that discontinuity flips direction as frequency changes from slightly below to slightly above an exact integer periodic-in-aperture frequency.  This is part of the effect of the default rectangular windowing of any finite length DFT or FFT.  Thus making the FFT phase result of a frequency sweep hard to interpret.

The way to remove this measurement at a circular discontinuity is to shift the waveform discontinuity away from the phase measurement point (FFT input sample 0) by doing a (circular) FFTShift before the FFT.  Then your FFT phase result will be referenced to the center of the original input vector (original waveform sample N/2), which will be continuous for continuous input functions or waveforms.

After the FFT following an fftshift, you can use the FFT frequency and phase results for the the center of the original input waveform to compute the phase at any other point (for instance the beginning, or sample 0 of the original input waveform.  Or the end, etc.).


There are many ways to skin this cat. I'll explain one way that I know is used in industry.

Create candidate frames algorithmically using the image features; e.g., the histogram, sharpness, presence of humans, faces, and anything else you can think of.
Serve different candidates to different users and see which ones are popular
You can stop once the popularity contest results stabilize. However, for bonus points you can learn the ranking function that yielded the list based on the features enumerated in part one.
Use the ranking function to select the best frame, or cull the list of candidates for subsequent tests.


Determining height of object from images

If you are only interested in the measurement of a single dimension (such as the height or the width or radius) of an object then you do not need multiple images. A single image would have worked if you could determine some scale factors which would depend on the distance of the camera to the object and orientation angles (which would be defined in a scaling matrices) etc.

If the measured dimension was inclined in 3D wrt the camera, then you should also apply projective geometry corrections which would also be defined in respective matrices.

Two images, on the other hand, would let you create the complete 3D environment through which any measurement is easily accomplished by their respective 3D point coordinates inside the world coordinate system. However it's considerably more complex and sensitive to create a precise definition of 3D points from just two images due to several reasons. You need very strict control on the imaging setup.


bit reversed order of FFT matrix for channel estimation

So, writing your $F^b$ and $F^{'b}$ as permutations $p, p^{-1}$ of the DFT matrix $F$ , IDFT $F'$, and the cyclic prefixing and removal also as Matrix operation $P$, $R$ (which boils down to cyclically extended or truncated identity matrices), and the circular convolution with the channel impulse response as circulant matrix $C$,
\begin{align}
x&=F^{'b}X=p F' X\\
x_{cp} &= PpF' X\\
y_{cp} &= y*x_{cp} \\
y&=R(y*x_{cp})= CpF'X\\
Y&=F^By \\
&= p'FCpF'X & F, F' \text{ are symmetric matrices!}\\
&=p'FCF'pX\\
&=p'\underbrace{FCF'}_{\text{IDFT-Channel-DFT}}\underbrace{pX}_{\text{permuted symbols}}
\end{align}
So, all that's really happening here mathematically is that you permute the symbols before transmission, then de-permute them at the receiver. In between, you get normal OFDM. There's no advantage or disadvantage to the permutation w.r.t. to channel diagonalizability, but we've arrived here in your previous question already.
Now, for channel estimation, you'll want to insert pilots such that they make sense for the central $FCF'$ part – in other words, you need to think about where you'd want your normal OFDM system to have these pilots, and then shift these positions according to your $p$.
Of course, the channel you estimate left of your final reordering isn't the same as the channel you estimate with OFDM, it's a permutation. As shown above, it's however "just" a permutation, there's nothing in the maths here that changes anything about the channels properties at all, so if your results aren't the same as for straightforward OFDM, then you've got an implementation bug somewhere.


Perform non-power-of-two FFT using ARM CMSIS library

Here is some code that calculates a 1280 tap FFT based on 5 256 tap FFTs

%% 1280 FFT based on five FFTs of 256 each
n = 1280;
% Create a piece of noise
x = randn(n,1);
% calculate FFT using MATLAB native fft() function.
% We'll use this as a reference to prove it works
fx = fft(x);

% Break down into five signals of 256 points each, interleaved
p = x(1:5:end);
q = x(2:5:end);
r = x(3:5:end);
s = x(4:5:end);
t = x(5:5:end);
% FFT each of those. This is a 256 power-of-two standard FFT
fp = fft(p);
fq = fft(q);
fr = fft(r);
fs = fft(s);
ft = fft(t);
% Do five times periodic extension (just repeat it 5 times)
fp5 = repmat(fp,5,1);
fq5 = repmat(fq,5,1);
fr5 = repmat(fr,5,1);
fs5 = repmat(fs,5,1);
ft5 = repmat(ft,5,1);
% calculate the 1280 twiddle factors
k5 = (0:n-1)';
W5 = exp(-1i*2*pi*k5/n);
% assemble the result
fy5 = fp5 + W5.*fq5 + W5.^2.*fr5 + W5.^3.*fs5 + W5.^4.*ft5;
% calculate the error
ferror = fy5-fx;
fprintf('Error = %6.2f dB\n',10*log10(sum(ferror.*conj(ferror))./sum(fx.*conj(fx))));


Image transformation where a curve becomes a straight line

I would recommend a grid transform. Each point of the grid has a $(dx,dy)$ offset, and a point in the middle of a square is interpolated according to the control points. Specifically in your case, you can keep most of the points in place, except those which are near the curve you want to straighten.


Blank as in all pixels are white or an "empty" image?
For the first case, you can simply check if all the pixels are white by using something like if np.sum(img == 255) == (number of element in image, i.e., width*heigth).
(Assuming it is a 8-bit grayscale image).
For the second case, you can open the image using openCV, PIL or matplotlib and then check the dimensions of the result.


How is the energy conservation principle preserved in the case mentioned below?

Power is homogeneous to a second degree quantity. When the data is doubled, the power is quadrupled.

The energy here (classicaly) is defined in an Hilbert space as a squared $L_2$ norm, deriving from some inner product $\langle\cdot,\cdot\rangle$, assuming a real Hilbert space,  $\|x\|^2 = \langle x,x\rangle$. Hence:  $$\|x+y\|^2 = \langle x+y,x+y\rangle = \|x\|^2 + 2\langle x,y\rangle+  \|y\|^2 \,.$$

Now set $y=x$, you gracefully get $\|x+x\|^2 = 4\|x\|^2$.


Why do eigenvectors specify the direction of fastest and slowest changes in the image

The eigenvector corresponding to the largest eigenvalue of the autocorrelation matrix indicates the direction of fastest change, while the eigenvector of the smallest eigenvalue indicates the direction of slowest change. Nothing strange here.

What's confusing is the Uncertainty Ellipse (Figure 4.6, page 213). The direction of slowest change is longer because it's scaled by the square root of the inverse eigenvalue of A. Really what it represents is the direction of the highest uncertainty.

Uncertainty and speed of change are different sides of the same coin:

high uncertainty <--> slow change
low uncertainty <--> fast change

The purpose of all this is to help us find good match points in an image. An example of a good match point is shown in Figure 4.5b. The reason it's "good" is because it has a well defined optimum. This occurs when the image has fast change in more than one direction (i.e. both eigenvalues of the autocorrelation matrix are large).

Figure 4.5c shows a match point that doesn't have a well defined optimum (it's more like a  ridge). The fact that we cannot precisely point to the optimum means that there is uncertainty in the direction of the ridge. The autocorrelation function of this patch will have one large eigenvalue and one small eigenvalue. The large eigenvector is in the direction of fastest change (a direction perpendicular to the ridge), while the small eigenvector is in the direction of the ridge.


How to interpret the sign of a cepstral coefficient?

For the zero-th cepstral coefficient, the sign is probably always positive since it represents the average energy of the entire spectrum.

For higher-order coefficients, seeing that the question was first posted in 2017, does the OP maybe already have an answer? Would be interested to learn.


FFT analysis of binned data of different length

Let's say you have several data sets $x_i$ and $N_i$ is the length of the $i$th data set.

First you're right that the amplitude of the FFT output scales with its length. To normalize all data sets they have to be divided by $N_i$. (Note that depending on the implementation of the FFT there might be a scaling involved already, check the Doc or see this question)

To combine all data sets you should append $(\max_i N_i) - N_i$ zeros to every data set before taking the FFT. The frequency axes are then aligned correctly and you can calculate the average of all FFT outputs.

Btw, I would plot the absolute value of the FFT output, not the real part. And sometimes I find it useful to "smooth" the spectrum by averaging every $M$ neighbouring FFT bins. Last, I have assumed that the data of every set has been acquired at equidistant time steps, otherwise the FFT isn't applicable at all.


Multiplication of two time domain sequences corresponds to cyclic (circular) convolution of the DFTs of the two sequences, not to linear convolution as implemented in your code. So if $x[n]$ and $y[n]$ are two length $N$ time domain sequences, and $X[k]$ and $Y[k]$ are their respective DFTs, the following correspondence holds:

$$\text{DFT}_N\left\{x[n]\cdot y[n]\right\}[k]= \frac{1}{N}\sum_{m=0}^{N-1}X[m]Y[(k-m)_N],\quad 0\le k \le N-1$$

where

$$(k-m)_N=k-m+lN$$

with the integer $l$ chosen such that the result is in the range $[0,N-1]$.

Also take a look at this answer to a related question.

EDIT: You also can't use Y1 in the convolution because it is not the DFT but the magnitude of the DFT. You should circularly convolve fft(y) with fft(w).


To find the impulse repsonse using the difference equation

There is no real shortcut for computing the impulse response. Partial fractions is the standard way to do it. However, in the case of a second-order transfer function as in your example, the result can be found in tables of $\mathcal{Z}$-transform pairs. E.g., if you use the last two correspondence of this table, you can pretty quickly write down the result.

By the way, your result can't be correct (at least not for negative $n$), because you know that the system is causal, so its impulse response must be zero for $n<0$.


How to represent impulse function in 2D?

So dealing with generalized functions like the Dirac delta requires some care, and when dealing with N-dimensional versions you need to be very explicit with your notation to keep things straight.

I'll denote the 2 dimensional delta function in polar coordinates at the origin as ${}^2\delta(r, \theta) = {}^2\delta(r)$, since for the special case of the origin, $\theta$ doesn't matter.

For this derivation, ${}^2\delta(r)$ represents the following limiting sequence of functions (that are asymmetrical about 0):

$${}^2\delta(r) = \lim_{\epsilon \rightarrow 0^+} \dfrac{1}{\pi\epsilon^2} \quad 0 < r < \epsilon$$

Using this limiting sequence of functions, the 2-D delta function at the origin in polar coordinates has the following property of "integrating to 1" for the 2-D integration:

$$\int_0^{2\pi} \int_0^{\infty} {}^2\delta(r) \space r \space \mathrm{d}r \space \mathrm{d}\theta = 1$$

To separate the $r$ and $\theta$ portions of ${}^2\delta(r)$, so we can use a 1D Dirac Delta function, I'll use the following limiting sequence of functions (that are asymmetric about 0) for the 1D Dirac Delta function:

$${}^1\delta_a(r) = \lim_{\epsilon \rightarrow 0^+} \dfrac{1}{\epsilon} \quad 0 < r < \epsilon$$

Using this limiting sequence of functions, the 1-D delta function at the origin has the following property of "integrating to 1" for the 1-D integration:
$$\int_0^{\infty} {}^1\delta_a(r) \space \mathrm{d}r = 1$$

We can then equate the two above integrals to get the relationship between  $ {}^2\delta(r)$ and ${}^1\delta_a(r)$

$$\begin{align}\\
\int_0^{\infty} {}^1\delta_a(r) \space \mathrm{d}r &= 1\\
&= \int_0^{2\pi} \int_0^{\infty} {}^2\delta(r) \space r \space \mathrm{d}r \space \mathrm{d}\theta\\
&= \int_0^{\infty} {}^2\delta(r) \space 2\pi r \space \mathrm{d}r\\
\end{align}$$

so we have

$${}^2\delta(r) = \dfrac{{}^1\delta_a(r)}{2\pi r}$$

This differs from your desired result by a factor of $2$, because I did not allow the 1 dimensional delta function to be symmertic around $0$, which wouldn't make sense for a polar origin where $r$ cannot be less than $0$ and $\theta$ ranges in $[0, 2\pi)$.

If one allows $r$ to go negative and $\theta$ to range in $[0, \pi)$, then one can show for a 1 dimensional delta function that is symmetric around $0$
$${}^2\delta(r) = \dfrac{{}^1\delta(r)}{\pi |r|}$$

(note the absolute value bars!)

This is a good example of why delta functions require careful definitions and considerations when making statements.

To show the projections of a 2-D delta function are themselves delta functions, start with a suitable limiting sequence of functions for ${}^2\delta(x,y)$, like an infinitely thin and tall rectangle function both in x and in y, such that

$$\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} {}^2\delta(x,y) \space \mathrm{d}x \space \mathrm{d}y = 1$$

and the answer should fall out from the limiting sequence of functions you use.


In some physico-chemistry methods, baselines can be modelled by closed-form lineshapes. This is not typically the case for chromatography, and it can be even worse depending on the type of chromatography (gaz, liquid).

So you have to rely on some morphological properties that you believe valid for your type of signal. You can model the background by a given polynomial, assume derivability of the peaks, etc. And then fit, or optimize a cost function.

In a recent work, we have modeled gaussian peaks as sparse with sparse derivative. We add that peaks are typically positive. Then we suppose that the baseline is smooth enough to remain after the  use of a low-pass filter. Finally, we consider Gaussian noise, yet it works with Poisson noise as well.

The BEADS (Baseline Estimation And Denoising w/ Sparsity) method is described, and compared with two other methods in  "Chromatogram baseline estimation and denoising using sparsity (BEADS)", by Xiaoran Ning, Ivan W. Selesnick, Laurent Duval, in Chemometrics and Intelligent Laboratory Systems, December 2014, http://dx.doi.org/10.1016/j.chemolab.2014.09.014.

We have made the associated Matlab toolbox available in Matlab.  Contributed versions in C++, R, and Python are made available through the companion BEADS page.


How to Save FFT output for fast read-in?

That doesn't seem worth the bother.

FFT of a 7 second long wave file on my windows 10 laptop using Matlab takes about 5 milliseconds. So unless you have a particularly slow setup, you will not save a significant amount of time by storing the FFTs instead of just reading the raw wave file and doing the processing again each time.

That doesn't seem to justify the extra development effort for a non-standard file format.


Output of system defined by differential equation

What you want to do is take the Laplace transform of both sides of S.  The first and second derivative terms will be replaced by expressions of $Y(s)$ and the conditions at $t=0-$, so you can solve the equation for $Y(s)$ and then do the inverse Laplace transform, or more commonly, use partial fraction expansion to go back to the time domain.

See numbers 35 and 36 in the table of Laplace transforms here.


PSOLA can mean different things to different people.

the "Formant Preserving Pitch Shifting" that is done in the time domain requires a really good, tracking, pitch detector.  knowing the pitch (which changes) at every arbitrary time means you know the period length.

so this Formant Preserving Pitch Shifting in the time domain, identifies the onsets of every period of the input waveform, when it is voiced.  each period is isolated with a window and now becomes a "grain" or (my poorly-considered name for it) a "wavelet".  then you launch these grains or wavelets out at the rate of the new fundamental frequency.  but, to preserve the formants, you neither stretch (if down-shifting) nor scrunch (if up-shifting) the grains or wavelets.  you overlap and add them.

the rate of firing out those grains determines the output pitch.  the degree of stretching or scrunching the grains along the time axis determines the formant shift.


Calculating the true values of biquad filter coefficients from normalized form

Note that a biquad has $5$ degrees of freedom (not $6$), because $a_0$ can always be chosen as $a_0=1$ without loss of generality:

$$\begin{align}H(z)&=\frac{b_0+b_1z^{-1}+b_2z^{-2}}{1+a_1z^{-1}+a_2z^{-2}}\\&=b_0\cdot\frac{1+\hat{b}_1z^{-1}+\hat{b}_2z^{-2}}{1+a_1z^{-1}+a_2z^{-2}},\qquad \hat{b}_1=\frac{b_1}{b_0},\;\;\hat{b}_2=\frac{b_2}{b_0}\tag{1}\end{align}$$

Eq. $(1)$ shows that $b_0$ is just an overall gain (for fixed values of $\hat{b}_1$ and $\hat{b}_2$), as you've noted yourself.

If you don't know the formula that produced the filter coefficients, there is generally no way to know the intended gain. As mentioned in Dan's answer, the gain is often chosen such that quantization effects are minimized. It can often be compensated for after all computations have been performed.

In some cases, however, it is possible to estimate the intended gain. E.g., if the biquad is a low pass filter, it is reasonable to assume that the original frequency response has a value of $1$ at DC. This means that

$$b_0\frac{1+\frac{b_1}{b_0}+\frac{b_2}{b_0}}{1+a_1+a_2}=1\tag{2}$$

Similarly, for a high pass filter a common scaling makes sure that the frequency response equals $1$ at Nyquist, which implies

$$b_0\frac{1-\frac{b_1}{b_0}+\frac{b_2}{b_0}}{1-a_1+a_2}=1\tag{3}$$

You can come up with similar guesses for the intended scaling for other standard filter types, like band pass filters (unity gain at the center frequency) and band stop filters (unity gain at DC and/or Nyquist).


Missing known sinusoid in FFT of sensor data

If your sampling frequency is 5Hz (Which I don't think is!) and the signal contains a known frequency of 5Hz. So, the sampled signal will have aliasing for sure.

You are correctly mapping your positive FFT spectrum from 0 to 2.5Hz, because digital frequency of sampled signal can only map from $[-f_s/2, f_s/2]$ which is $[-2.5hz, 2.5hz]$.

The question is where do you want to see your 5 Hz in this. You can only see a high DC which is actually the negative folded image of the 5Hz tone at $f_s$. And, you can actually see a peak very close to DC.

  Based on the Plots you have shown :

The plot of time domain data and assuming you have injected a 5Hz tone, I think the sampling frequency is close to 500Hz and not 5Hz. Simply, because you have ~2.5 cycles of 5Hz sinusoidal in 250 samples. 5Hz sinusoidal wave will take 1 second for 5 cycles, hence ~2.5 cycles of such wave will take 0.5 second. 250 samples/0.5 second means 500 samples/second. Hence, your sampling frequency is 500 Hz.

And, if you are taking 250 length FFT, frequency resolution will be of :
$$\frac{2\pi}{250} = \frac{2}{250}.\frac{f_s}{2} = 2Hz$$
Hence, your 5Hz tone will be seen as a peak between $2^{nd}$ and $3^{rd}$ coefficient of the 250-Point DFT, which you are able to see.

MATLAB Simulation showing your situation:

fd = 5;             % data tone of 0.25Hz
fs = 500;           % sampling frequency of 5Hz
Ts = 1/fs;
n_samps = 250;      % Number of samples
time = n_samps/fs;  % duration of data collection (0.5 secs)
n = linspace(0,time-Ts,n_samps); % sample points in the duration [0sec, 0.5sec)

% Carrier Frequency fc = 5Hz + a slow linear ramp as it seems from the
% time-domain plot you have attached.
signal = (n/250 + sin(2*pi*fd*n));

% AWGN Noise added
signal = awgn(signal, 15);
figure(1)
plot(signal)

Y = fft(signal);
figure(2)
plot(5/n_samps * (0:1:n_samps/2-1), (abs(Y(1:length(Y)/2))));

  The other small peaks are most probably being generated due the
  harmonics which are present as a consequence of clipping in sampled
  data due to fixed bit-width ADC. If you lower your amplitude a little
  to accommodate your sinusoidal wave within the dynamic range of ADC,
  these small peaks might go away.

Adding analysis of the time-domain samples you have pasted :

Assumption : Sampling Frequency = 500Hz

I have plotted FFT of the samples centered around 0 Hz(DC). And, based on number of samples (256) and assumed $f_s = 500Hz$, frequency spacing between 2 FFT samples is $\frac{f_s}{n_samps} = \frac{500}{256} = 1.9531$. Since 5Hz is not an integral multiple of this spacing, therefore, you will see the tone energy spread around into 1.9531*2Hz and 1.9531*3Hz i.e. 2nd and 3rd FFT coefficient.

You are also seeing a strong peak at DC (0Hz) because your sinusoid of 5Hz is having an amplitude shift of a non-zero value.

You can see that your tone is not missing at all. It is present right where is should be.

MATLAB Code to plot the signal and it's DFT with correct $f_spacing$ is as follows:

plot(500/256*(-128:127), abs(fftshift(fft(x))))


How to do prediction using frequency domain data?

An important NOTE: Since you are talking about frequency
  domain, it is implied that entire DFT spectrum is available and hence
  estimation is used for smoothing rather than future prediction.

If the signal is stationary, you can apply wiener filter and the model produced is an FIR filter; in this case, the signal estimation in the time domain will be identical to that of frequency domain.

From wiki:
Wiener's main accomplishment was solving the case where the causality requirement is in effect, and in an appendix of Wiener's book Levinson gave the FIR solution.

Removing noise using wiener filter using deconvolution is called Wiener deconvolution. This works in frequency domain. And is quite well used in Image deconvolution.

I don't know if there is a formulation possible for Kalman filter to be used for given Frequency domain data (assuming DFT) because usual implementations are actually iterative sample by sample.  But kalman smoothing approaches probably can do similar thing.


The signals plotted on the figure are correct.

As your provided link already shows: $$y_d(t) = \mathcal{T}\{x_d(t)\} = \mathcal{T}\{x(t-d)\} = \mathcal{T}\{x(t-2)\} = x(2t-2) \not= y(t-2) $$ where $$ ~ y(t)=\mathcal{T}\{x(t)\} = x(2t) \longrightarrow y(t-2) = x(2(t-2)) = x(2t-4) $$

Looking at the progression on the provided figure, $x_1(t)$ is a pulse of domain $~ [-2,2] ~$ then $y_1(t)=\mathcal{T}\{x_1(t)\} = x_1(2t)$ will be a pulse of domain $~ [-1,1] ~$ which is obtained by compressing $x(t)$ by $2$ about the origin $~t=0~$

$x_2(t)$ is obtained from $x_1(t)$, delaying it by $2$ units to the right, hence $x_2(t)$ has a domain of $~[0,4]$.
Therefore $y_2(t)=\mathcal{T}\{x_2(t)\} = x_2(2t)$ will be a pulse of domain $~ [0,2] ~$ which is obtained by compressing $x_2(t)$ by $2$ about the origin $~t=0~$

On the other hand shifting the pulse of $y_1(t)$ with the domain $[-1,1]$ by $2$ units to the right creates a new pulse with domain $[1,3]$ which is not equivalent to $[0,2]$; the domain of $y_2(t)$, from which we deduce, also graphically, that $y_1(t-2) \not = y_2(t)$


Measuring similarity of music, using lossy compression

Is there a good way to use them to measure how similar the two clips
  are?

This is a very sane intuition but something I have rarely seen in the literature - outside the peripheral idea of doing feature extraction from compressed streams (and the motivation here is just the reduced computational burden).

I think the major reason is that audio codecs do not store a perceptual representation of the audio data. Instead, they store a very coding-centric representation of audio - entropy-coded, quantized, coefficients of a transform - chosen so that the distortion between the original and the encoded signal falls below the threshold of perception. In a sense, looking at a compressed audio stream does not tell you anything about how a human would perceive the audio signal. It just tells you that anything that would fall below the perception thresholds has been taken away. The codec "carves" around the meaningful audio signal, but never touches it.

  where L(x) is the length of the compression of sound clip x, and C1||C2 is the result of concatenating the two clips?

This makes sense on strings - I think this idea is found in some papers by Cilibrasi or Vitanyi - but this is not how audio codecs work. The audio codecs you mention are designed to process streaming audio with low latency, so they only "see" a short window of the signal at a time and do not attempt to eliminate long term redundancy. Just like JPEG for images. Copying 10 instances of the same note one after the other will result in a file 10x as large. Audio codecs have a fixed bit-rate, or an adaptive bit-rate that is adapted according to the content of the narrow window of signal that the codec sees.

The redundancy eliminated by an audio codec is very short-term.

Your idea would require a much more sophisticated "objet oriented" codec that attempts to decompose the audio into a hierarchy of objects, such as notes. This idea is still in its infancy in the academic world.

Another way to look at it is to check if the decoding process bears similarity with a generative model of the signal. This is the case for some string compression algorithms, where the compression process recovers something akin to a generative grammar of the string (Sequitur, factor oracles...). This is the case for LPC-based speech codecs in which the LP coefficients embed knowledge about the articulation, and the residual knowledge about prosody and voicing/unvoicing - the compressed stream can be seen as data for a speech synthesizer. But this is not the case for general purpose audio codecs like the ones you mention - the decoder bears no similarities with the music production process.

  Or maybe find the highest bitrate such that F(C1) is close to F(C2)) via a simple metric (e.g., L2 norm applied to the FFT spectrum)

You might do the opposite - look at the bitrate under which the clips become indistinguishable, but ultimately it wouldn't be that much easier than comparing sinusoidal representations pruned by a perceptual function.


Why Fourier transform is not sufficient and we have to use Laplace transform?

The Laplace Transform is more representative of real systems that have a starting point, which is why the integral starts at 0, and also why the unit step function is generally talked about alongside the Laplace Transform. With the Laplace Transform, we can examine the transient and steady-state behavior of a system.

Using $e^{st}$ instead of $e^{iwt}$ allows us to examine different aspects of a physical system. The variable $s$ is complex, and if the real part was set to 0, it would reduce to a truncated Fourier Transform. The real part of $s$ is related to the amount of damping in the system. Also, with the Laplace Transform, a system's stability can be considered.

In short, Laplace is used to consider damping, stability, transient and steady-state behavior of a physical system (represented by a differential equation).


Linear convolution of discrete signals with defined lengths

It seems like you have already the correct answer, but try to visualize what's going on

First understand that signals of length $n_0$ are really infinite length, but have nonzero values at $n = 0$ and $n = n_0 - 1$. The values in between can be anything, but for the purposes of this problem take them to be nonzero as well.

Now perform the discrete convolution by literally shifting the length-5 signal and dot multiplying it with the length-7 signal. Your result will also be an infinite length signal with nonzero values only where the two signals overlap (when they dont overlap, you should find the convolution to be zero). In this case, there are 11 points where the signals actually overlap.

If some parts within the signal are zero, it is possible that you get fewer nonzero values in the result. However, in the max case where the full signal is nonzero you get this max, $11 = 7 + 5 - 1$ samples


Practical insightful time domain functions

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

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

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


Window function for signal with 3-5 sine waves inside

Windowing, in the context of DFT, is used to reduce the side-lobes or spectral leakage which happens because you work with finite length discrete signals.

Now, if is sum of sinusoids that are periodic in the observation time (that means that each sinusoid completes an integer amount of periods), then the DFT will be a series of delta function with no side-lobes.
However, in a more realistic scenario, the DFT of your signal should be effected by the spectral leakage. In that case, it is possible that some of the side-lobes will have greater magnitude than some of the sinusoids and if you try to select the 5 frequencies with the greatest magnitude, you will get the wrong results.
To sum it up simply: Yes, you should use some kind of windowing. It doesn't actually meter which window function you use. As suggested, Hamming window is a good choice.


Why is it necessary to average signals with filters?

An average is a summary of different values: if you have several signal processing exams, you know that both your answers and the grading person may incur some variability. On top of that, it is complicated to compare student, or the same student at different years, based on a series of grades.

So an average of grades both:

compacts different grades in one single number,
balances exam variability.

The latter is of first importance in signal processing: a signal sample may be disturbed by noise, sensor fluctuations, quantization, etc. While the signal is supposed to have some regularity, the disturbances are often chaotic, incoherent, from one sample to the other. If you average the distances alone, their amplitudes tend to diminish (not in all cases, but let us say that for now). So if close sample values are averaged, the ones from the signal don't change a lot (like a regular working student), while the disturbances reduce. Thus, the signal looks smoother, more pertinent. This is the basis for the linear filters: further, you can weight the samples differently.

The summarizing property is also important. You can average three samples, and keep the average only, and switch to the other sample triplet.  This reduces the quantity of information, and is the basis for subsampling, data compression, etc. JPEG uses a lot of filtering and subsampling.

Both combined are the root of filter banks, putting several filters in parallel with subsampling. Then, you have nonlinear filters, but this is a whole new story.


The idea of the algorithm is this: assume you have a length $N$ signal that is sparse in the frequency domain. This means that if you were to calculate its discrete Fourier transform, there would be a small number of outputs $k \ll N$ that are nonzero; the other $N-k$ are negligible. One way of getting at the $k$ outputs that you want is to use the FFT on the entire sequence, then select the $k$ nonzero values.

The sparse Fourier transform algorithm presented here is a technique for calculating those $k$ outputs with lower complexity than the FFT-based method. Essentially, because $N-k$ outputs are zero, you can save some effort by taking shortcuts inside the algorithm to not even generate those result values. While the FFT has a complexity of $O(n \log n)$, the sparse algorithm has a potentially-lower complexity of $O(k \log n)$ for the sparse-spectrum case.

For the more general case, where the spectrum is "kind of sparse" but there are more than $k$ nonzero values (e.g. for a number of tones embedded in noise), they present a variation of the algorithm that estimates the $k$ largest outputs, with a time complexity of $O(k \log n \log \frac{n}{k})$, which could also be less complex than the FFT.

According to one graph of their results (reproduced in the image below), the crossover point for improved performance with respect to FFTW (an optimized FFT library, made by some other guys at MIT) is around the point where only $\frac{1}{2^{11}}$-th to $\frac{1}{2^{10}}$-th of the output transform coefficients are nonzero. Also, in this presentation they indicate that the sparse algorithm provides better performance when $\frac{N}{k} \in [2000, 10^6]$.

These conditions do limit the applicability of the algorithm to cases where you know there are likely to be few significantly-large peaks in a signal's spectrum. One example that they cite on their Web site is that on average, 8-by-8 blocks of pixels often used in image and video compression  are almost 90% sparse in the frequency domain and thus could benefit from an algorithm that exploited that property. That level of sparsity doesn't seem to square with the application space for this particular algorithm, so it may just be an illustrative example.

I need to read through the literature a bit more to get a better feel for how practical such a technique is for use on real-world problems, but for certain classes of applications, it could be a fit.


I assume you mean that at each time $t$, the signal is a normally-distributed random variable.  This tells you the probability that the signal will take a value in any given range, but it does not tell you whether/how the signal values at different times are related to one another.  Typically, for a random signal one defines not only the probability distribution for the signal at each time, but also the autocorrelation function $E \{ x(t) x(u) \}$ or some other measure of how signal values at different times are related or independent.  This is required in order to determine how the signal is distributed in the frequency domain.


Each frame of your movie is a 2D plane indexed by $x$ and $y$.
When you consider several consecutive frames, you can stack them in order to have a 3D cube, where each slice of the cube is a frame.
This space-time volume is now indexed by $x,y$ and $t$.

You can manipulate this volume as any 3D volume, just replace i nyour mind $t$ by $z$ if you find it confusing.

Once you have this volume, if you change the normal axis used to define the 2D slices (remember it was the $t$ axis initially), you obtain the different slices that you want:

choosing the $x$ axis as normal, you obtain 2D $(y,t)$ slices
choosing the $y$ axis as normal, you obtain 2D $(x,t)$ slices

By the way, the classical term for these kind of slices is epipolar-plane image analysis if you are loking for more references.


When you use the channel matrix (it is probably a convolution matrix), you just need to multiply it by the signal. Note that the output from this operation will not be equal the output from the filter function. The filter function generates an output from the same size of the input signal, so if your signal has a length of N, just the first N samples obtained from the multiplication will be the same. See the code posted by @Harris.


linear-expansion filter for image (2D array)

Simply resize the binarized image to a larger image but do not use any high order interpolation in resizing and instead use 'nearest neighbor'interpolation technique. In MATLAB:

im2 = imresize(image, 10, 'nearest')

part of 10x expanded image:

and a raster of the image:


Given an image how should I detect this is a day picture or night picture?

This will work only if you're dealing with natural scenes (no artificial light sources).

I think there are two other possible ways:

1) looking for the percentage of dark pixels (define a threshold) that will tell you how many black areas are present in the image. Again, this will work only if we're not considering photos close to artificial light sources. Doing an average of the pixel value will lower the result if you have many lights in the scene, considering only dark areas could lead to better result (if there are enough in the scene).

2) Natural light have some difference in frequencies than artificial light, you could search which one is more present in the photo. I'm thinking at night pictures of cities for example, where lights are often in the yellow range. But i think this seems easier to say than implement, and looks more conditioned to the type of photo (natural panorama, city images, etc.).

P.S: i'm assuming the picture has been taken with the correct exposure time


In line with a previous similar question here are my suggestions:

There are so many nice books but I believe you should first have a look at the science of sound from Rossing for getting the most broad view on the subject.

Then you can look at the following books , each dealing with a particular dimension of the problem:

1- Acoustics_BERANEK

2- Elements of Acoustical Engineering_OLSON

3- Music, Physics and Engineering_OLSON

4- Art of Digital Audio_WATKINSON

5- Audio Engineering_SINCLAIR

6- Digital Audio Effects_ZÖLZER

7- Digital Audio Signal Processing_ZÖLZER

8- Fundamentals of Acoustic Signal Processing_TOHYAMA

9- Audio Signal Procesing and Coding_SPANIAS

10- Audio Programming Book_BOULANGER

11- PsychoAcoustics_ZWICKER

12- Acoustics and PsychoAcoustics_HOWARD

EDIT:
I'm sorry that I've forgotten to mention about books on speech processing, so here you can find a list:

1- Speech and Audio Signal Processing_GOLD

2- Speech and Hearing for Communication_FLETCHER

3- Discrete Time Processing of Speech Signals_PROAKIS

4- Fundamentals of Speech Recognition_RABINER

5- Introduction to Digital Speech Processing_RABINER


Drawing a Spectrogram off of FFT - what are the decibels relative to?

For audio files you want to use dBFS, dB relative to full scale:

https://en.wikipedia.org/wiki/DBFS


importance of using histogram equalization

The difference is that with your method, if you imagine the histogram, you're simply going to stretch it to span from 0 to 255, but its shape will be preserved. Histogram equalization not only stretches your histogram, but also tries to make it flat, so that you get an approximately even distribution of pixels of every shade of gray.

In terms of why one is better for certain applications than others, that's application specific.

Edit:

Here are some examples from Bruzed:

Contrast stretching (your method):

Histogram equalization:

As you can see, the bulk of the pixels in the original image was gray, represented by a large peak in the middle. When you do contrast stretching, the peak is still there, even though your darkest pixel is now black and you brightest pixel is now white. By contrast (pun intended), using histogram stretching you get a much flatter histogram response. This actually increases the overall contrast of the image.


Bandpass filter with very high sampling rate

Limited numerical precision makes your filter unstable.

In this case you can fix it by implementing the filter in second order sections using the zeros and the poles of the filter.


Phase difference measurement of a signal sampled with two different sampling frequencies

This is an untested time domain solution, but the math looks solid.

This will be impossible to implement unless you solve the reciever synchronization problem first.  That is either a hardware fix or a calibration operation.  Assume it is solved and your two signals are coming in as time aligned sequences.  Assume also your sampling rates (I don't like "sampling frequency" as a term) are high compared to the frequency of the incoming signal.  Also the incoming signal is assumed to be a pure real tone of a known frequency ($\omega$).

You will be searching for peak values on each sequence.  Use a parabolic fit with the center point being the max.  You will then need to find the signal values from the sequences at that time instance by interpolation.

Another problem you will have is normalizing the magnitude of one signal against the other.  This solution does that implicitly.

With that in mind, we can do the math with a continuous definition.

$$ x_1(t) = A_1 \cos( \omega t + \phi ) $$

$$
\begin{aligned}
x_2(t) &= A_2 \cos( \omega ( t - d ) + \phi )  \\
&= A_2 \left[ \cos( \omega t + \phi ) \cos( \omega d ) + \sin( \omega t + \phi ) \sin( \omega d )  \right] \\
\end{aligned}
$$

Divide the second signal by the first.

$$ \frac{x_2(t)}{x_1(t)} = \frac{A_2}{A_1} \left[ \cos( \omega d ) + \tan( \omega t + \phi ) \sin( \omega d ) \right] $$

At a peak of signal 1, $ \tan( \omega t + \phi ) = 0 $, since the $\sin$ does which is the derivative.

$$ (\frac{x_2}{x_1})_{peak1} = \frac{A_2}{A_1} \cos( \omega d ) $$

If your amplitudes are normalized, you have your answer from this.

By symmetry the reverse argument can be made.

$$ s = t - d $$

$$ x_2(s) = A_2 \cos( \omega s + \phi ) $$

$$
\begin{aligned}
x_1(s) &= A_1 \cos( \omega ( s + d ) + \phi ) \\
&= A_1 \left[ \cos( \omega s + \phi ) \cos( \omega d ) - \sin( \omega s + \phi ) \sin( \omega d )  \right] \\
\end{aligned}
$$

$$ \frac{x_1(s)}{x_2(s)} = \frac{A_1}{A_2} \left[ \cos( \omega d ) - \tan( \omega s + \phi ) \sin( \omega d ) \right] $$

$ \tan( \omega s + \phi ) $ will be zero at a peak of signal 2.

$$ (\frac{x_1}{x_2})_{peak2} = \frac{A_1}{A_2} \cos( \omega d ) $$

We can now combine the ratios from the two different peak locations:

$$ \cos( \omega d ) = \sqrt{ (\frac{x_2}{x_1})_{peak1} \cdot (\frac{x_1}{x_2})_{peak2} } = V $$

Since we are squaring, the sign of $V$ is lost.  It can be deduced from the individual parts.

$$ d = \frac{ \cos^{-1} (V) }{ \omega } $$

Note, the solution is not unique as the shift can be greater than one cycle length.

This should take a lot fewer calculations than an FFT approach.

Had a "duh moment."

If you can measure the peaks directly, then $d$ is simply their distance apart plus possibly a number of whole cycles.

You can calibrate your system if you can move your source straight ahead for a baseline shift measurement, ideally zero.  Then a measurement the same distance on either side will give you an idea of the angle of the first cycle overlap if there is one.

For a third solution, there is a simple frequency domain approach that will work well in this situation if your signal is too noisy for the time domain one.  Select a duration which is a whole number of cycles, say $k$, of the signal you are receiving.  It will be roughly $M$ number of samples at your one rate and $N$ on the second.  Calculate the $k$th DFT bin for each using two sets of basis vectors (two different DFTs, technically, but one bin each).  You should be able to find an interval so the $M$ and $N$ are pretty good fits.  Read the phase difference from the difference of the two bin angles.  Translate that to a time shift using the frequency.  This should still take a significantly fewer number of calculations than doing FFTs (DFTs) on each signal.


Concatenated Deletion & Sticky Channel Capacity

The identical and identically distributed (IID) binary deletion and sticky channels are both general repeat channels each characterized by a discrete probability distribution of repeating an input bit in the output $D$ times with $D \in 0, 1, 2, 3, 4, \dots$ with $D\!=\!0$ meaning that the input bit does not appear in the output of the channel at all. The question's probabilities for the deletion channel (first channel, hence the subscript $1$) are given by:

$$P(D_1) = \left[\tfrac{1}{2}, \tfrac{1}{2}, 0, 0, 0, \dots\right],\tag{1}$$

which is to say $P(D_1\!=\!0)=P(D_1\!=\!1)=\frac{1}{2},$ as given in the question.

An IID sticky channel with input duplication probability $\frac{1}{2},$ as given in the question, and of a variety that may at most double the number of bits, is also a general repeat channel, with $D_2\text{-}$times repeat probabilities:

$$P(D_2) = \left[0, \tfrac{1}{2}, \tfrac{1}{2}, 0, 0, \dots\right]\tag{2}$$

Generally, without hard-coding the above numbers, the concatenation (cascade, in series) of the above two channels into a single channel is also a general repeat channel with $D_1D_2\text{-}$times repeat probabilities:

$$P(D_1D_2) = \left[P(D_1\!=\!0) + P(D_2\!=\!0)\\- P(D_1\!=\!0)P(D_2\!=\!0),\\
P(D_1\!=\!1)P(D_2\!=\!1),\\
P(D_1\!=\!2)P(D_2\!=\!1) + P(D_1\!=\!1)P(D_2\!=\!2),\\
P(D_1\!=\!3)P(D_2\!=\!1) + P(D_1\!=\!1)P(D_2\!=\!3),\\
P(D_1\!=\!4)P(D_2\!=\!1)\!+\!P(D_1\!=\!1)P(D_2\!=\!4)\!\\+\!P(D_1\!=\!2)P(D_2\!=\!2),\\
\dots\right]\tag{3}$$

These are probabilities for the possible values of the product of independent variables $D_1$ and $D_2.$ For example the last entry before the ellipsis collects probabilities for the three mutually exclusive events $4\times1=$ $1\times4=$ $2\times2$ $=4.$
Note that $P(D_1D_2) = P(D_2D_1),$ so it does not matter which channel comes first. Using the values of $P(D_1)$ and $P(D_2)$ from Eqs. 1 & 2, we get:

$$P(D_1D_2) = \left[\tfrac{1}{2}, \tfrac{1}{4}, \tfrac{1}{4}, 0, 0, \dots\right]\tag{4}$$

This characterization of the concatenated channel appears more complex than the same for the deletion and sticky channels. For each of the latter, finding the capacity appears to be an open research problem (now at the end of 2018). So I doubt the capacity of the concatenated channel is known either. However, it may be possible to find upper and/or lower bounds. For the capacity of the deletion channel, there seems to be a lower bound $1/18,$ see: Capacity of a Binary Deletion Channel.

Your question says that there should be no error correction. As far as I understand, the channel capacity does not depend on error correction, but if a working code can be demonstrated, then it gives a lower bound for channel capacity.

Mahdi Cheraghchi calls the general repeat channels $\mathcal{D}\text{-}$channels, where $\mathcal{D}$ would be the distribution of $D=D_1D_2$ for the concatenated channel. Perhaps you will find something of interest in their work (Cheraghchi 2018):

  A $\mathcal{D}\text{-}$repeat channel, given a bit,
  draws an independent sample $D ≥ 0$ from $\mathcal{D}$
  and outputs $D$ copies of the received bit.

  Our framework can be applied to obtain capacity upper
  bounds for any repetition distribution (the deletion and
  Poisson-repeat channels corresponding to the special cases
  of Bernoulli and Poisson distributions).

References

Mahdi Cheraghchi. 2018. Capacity upper bounds for deletion-type channels. In Proceedings of the 50th Annual ACM SIGACT Symposium on Theory of Computing (STOC 2018). ACM, New D_2ork, ND_2, USA, 493-506. DOI: https://doi.org/10.1145/3188745.3188768 (arD_1iv:1711.01630 [cs.IT])


HMM-based speech recognition/synthesis: state representation of monophones and triphones

My current understanding is that each of the phones in a triphone corresponds to each of the states.

That's not correct.

Each phone, regardless of the context dependency, has three states: "start", "in-between", and "end". There are always transitions from the previous state to the next one and to itself. This defines the structure of HMM. Now context-dependency means that for each triphone we know its left and right context (context-independent phone). Let's take a look at what phone AX could look like in both CI and CD HMMs for the word "one [W AX N]".

When context-independent phones are used AX is begin split it into three states - _AX_B_, _AX_I, _AX_E - to construct the corresponding HMM:

In case of context-dependent phones AX actually doesn't exist any more. Instead we have many new phones that store information of phones to the left of it and to the right of it. Thus in the given word the neighbours of AX are W(left) and N(right), so HMM for phone W-AX-N is constructed:


The meaning of the sign in a phase difference

I believe the terms “leading” and “lagging” can be misunderstood, so I don’t recommend using that terminology.

Due to the traditional definition of the discrete Fourier transform (DFT), a measured DFT sample phase of zero degrees means the input sinusoid was a cosine wave with an initial phase of zero degrees at time t = 0. So when you measure a DFT sample’s phase to be –90 degrees (your ‘A’ signal) that means the input sinusoid was a cosine wave with an initial phase of –90 degrees which is equal to a sine wave with an initial phase of zero degrees at time t = 0.

Your signal ‘B’ phase of –120 degrees means that input signal was a sine wave with an initial phase of -30 degrees at time t = 0.