Analysis is separate from processing.  Windowing helps remove certain artifacts during spectrum analysis, but since it is numerically lossy, isn't as useful when actually processing the audio for output.  Instead, sufficient zero-padding (without any non-rectangular windowing) and the use of overlap-add or overlap-save is used for any fast convolution or frequency domain filtering.  This prevents (actually just zeros) circular convolution artifacts.

Convolution reverb is already answered elsewhere, but in general can be done by breaking up the reverb impulse response into a whole series of overlap-save processing steps and summations (e.g. a 256k long reverb might be done using a series of 1,000 256+padding length FFT convolutions, or a mix of shorter and longer ones, with a future use buffer to hold most of those results till needed.)


If I understand correctly, the OP is comparing cascading 5 time domain signals one after the other as in $[f_1, f_2, f_3, f_4, f_5]$ and taking the DFT to summing longer samples of the signals and then taking the DFT. The following is what we would expect to see in this case.

The DFT of the sum of functions is equal to the sum of their DFT's:

$$\sum_{n=0}^{N-1}f_1(n)W_N^{nk}+\sum_{n=0}^{N-1}f_2(n)W_N^{nk} = \sum_{n=0}^{N-1}(f_1(n)+f_2(n))W_N^{nk}$$

Consider first all sequences of the same length and cascading versus summing the sequences:

Cascading $R$ sequences of equal length $M$ is identical to zero-padding each of the sequences out to $N = RM$, each with  starting offset of $rM$ where $r$ is from $0$ to $R$.

Example with two sequences with $M=4$:

$f_1 = [1, 2, 3, 4]$

$f_2 = [6, 8, 9, 1]$

$\text{DFT}\{[1, 2, 3, 4, 6, 8, 9, 1]\}$

$ =  \text{DFT}\{[1, 2, 3, 4, 0, 0, 0, 0]\} + \text{DFT}\{0,0,0,0, 6, 8, 9, 1]\}$

Zero-padding a sequence prior to taking the DFT does not change the original DFT values, but inserts new interpolated values in between the samples in frequency. Such that

$$F_p(nk) = F(k) \tag{1}\label{1}$$

Where $F$ is the DFT of a sequence $f$ and $F_p$ is the DFT of $f$ with zero-padding. This is proven below along with the fact that equation $\ref{1}$ holds true for all offsets $rM$.

So in the example with $f_1$ aand $f_2$ above, $$F[k] = \text{DFT}\{[1, 2, 3, 4, 6, 8, 9, 1]\} $$

$$F[nk] = \text{DFT}\{[1, 2, 3, 4]\} + \text{DFT}\{[6, 8, 9, 1]\}$$

Similarly, the DFT of the cascade of the five sequences would have the same values as the DFT of the sum of the five sequences with four interpolated values in between each of the DFT samples!

Since in this case the OP didn't simply sum the sequences, but rather summed longer sequences out to the same length as the cascaded sequence, then here the DFT of the longer sequences would have the higher frequency resolution and scaled 5 times larger given by the longer data set, so the values in between that were previously interpolated would now be dictated by the longer time samples as new DFT values. However every 5th sample in both sets should match exactly other than the scale of 5 for the same reason outlined above, with more samples of actual frequency data in between!

Proof of relationship used above in equation $\ref{1}$

Given sequence $f(n)$ where $f(n) = 0 \text{ for } n \notin [0, M-1]$

The DFT of $f(n)$ is:
$$F(k)=\sum_{m=0}^{M-1} f(m)W_M^{mk}, \space\space\space\space k = 0,...M-1\tag{2} \label{2}$$

Where $W_M^{n}$ are the $n$ roots of unity given by $e^{-j2\pi n/N}$

The zero padded DFT of $f(n)$  out to length $N=RM$ with $R\ge 0$ with any offset $rM$ with $r \in [0, R]$ is:

$$F_P(k) = \sum_{n=0}^{N-1} f(n-rM)W_{N}^{nk}, \space\space\space\space k = 0,...N-1\tag{3} \label{3}$$

Using the DFT shift theorem, equation $\ref{3}$ becomes:

$$F_P(k) = \sum_{n=0}^{RM-1} f(n-rM)W_{RM}^{nk}=W_{RM}^{rMk}\sum_{n=0}^{RM-1} f(n)W_{RM}^{nk} $$
$$= W_{R}^{rk}\sum_{n=0}^{RM-1} f(n)W_{RM}^{nk}$$

Since $f(n) = 0$ for $n>M-1$

$$F_P(k) = W_{R}^{rk}\sum_{n=0}^{M-1} f(n)W_{RM}^{nk}\tag{4} \label{4}$$

$$F_P(Rk) = W_{R}^{rRk}\sum_{n=0}^{M-1} f(n)W_{RM}^{nRk} = \sum_{n=0}^{M-1} f(n)W_{M}^{nk} = F(k)$$

Thus we see how given any sequence $f(n)$ of length $M$ with DFT $F(k)$, zero padded out to length $N=RM$ for any positive integer $R$ with DFT $F_P(k)$, the following equality will hold for any rotational offset (starting index) of $rM$ for any integer $r$ from $0$ to $R-1$:

$$F(k) = F_p(Rk) \tag{5} \label {5}$$


A practical approach with centering, filtering, derivative, etc. may provide some results. Taking a step back on associated properties, for instance in the Fourier domain, may give you different insights, since symmetry is not well-defined per se (local or global, radial or axis, etc.). I am thinking specifically of the phase-based approach in Peter Kovesi, Symmetry and Asymmetry from Local Phase, 1997:

  Symmetry is an important mechanism by which we identify the structure
  of objects. Man-made objects, plants and animals are usually highly
  recognizable from the symmetry, or partial symmetries that they often
  exhibit. Two difficulties found in most symmetry detection algorithms
  are firstly, that they usually require objects to be segmented prior
  to any symmetry analysis, and secondly, that they do not provide any
  absolute measure of the degree of symmetry at any point in an image.
  This paper presents a new measure of symmetry that is based on the
  analysis of local frequency information. It is shown that points of
  symmetry and asymmetry give rise to easily recognized patterns of
  local phase. This phase information can be used to construct a
  contrast invariant measure of symmetry that does not require any prior
  recognition or segmentation of objects.

Here are examples for testing local symmetry:

And some Matlab code, for instance phasesym.m in Peter Kovesi Matlab functions.


Why is it easier to design an IIR filter in the continuous time domain than in the discrete time domain?

The term “easy” can mean simply to need less effort but can include cases where an accumulation of prior efforts can be leveraged.

Analog filters have a long  history and someone in 1970 who was trained in analog design might consider modifying what they know to build a digital filter as “easier" than using an optimization routine written in Fortran running on their IBM mainframe.


Explanation of FIR and IIR Filters with Simple Words

Filters modify the frequency content of a signal.

The terms Finite Impulse Response (FIR) filter and Infinite Impulse Response (IIR) filter refer to Linear Time Invariant (LTI) filters.

LTI filters can change the amplitude and phase of specific frequencies in a signal, but they cannot shift any frequency content to other frequencies.

LTI filters are usually used to remove undesired frequencies from an input signal.

That's it.

As an example, your FM radio receives a signal that is the entire FM broadcast band, from 88 MHz to 108 MHz in the electromagnetic frequency spectrum, at its antenna.  In a simplified view of an FM radio receiver, an LTI filter is used to filter out everything but the 200 kHz around 99.5 MHz, so that the radio can demodulate the audio from the station at 99.5 MHz without the frequencies from all the other stations interfering with that process.  (In actuality FM radio receivers are more complicated, using a cascade of LTI filters with some non-linear "mixing" processes, but we're going for the concept here.)

To describe some of the math, at a high level:

When FIR and IIR filters are applied to signal in time domain, their application is equivalent to a mathematical operation called "Convolution".

The mathematical operation of convolution of a signal with an FIR or IIR filter, is equivalent to time-reversing the filter response, shifting it relative to the signal, performing pointwise multiplication, then summing up the values; then shifting the time-reversed filter response over slightly relative to the signal, preforming pointwise multiplication, then summing up the values; then ... keep repeating that for all of time (or until the input signal disappears to 0).


As I've mentioned in a comment, the poles are not "ordered" in any way. You can come up with different versions of the formula for the pole locations, which will give you the same poles in different orders.

But anyway, looking at the formulae in your question, the first thing to notice is that there should be a negative sign in the formula for the real parts $\sigma_k$. Second, the values that you got must be the imaginary parts $\omega_k$ (and not $\sigma_k$ as written in your question). If I evaluate the formula for $\omega_k$ given in your question for $k=1,2,3$ I get the following values (for $n=3$ and $1$ dB pass band ripple):
$$\begin{align}\omega_1&=0.966i\\
\omega_2&=0\\
\omega_3&=-0.966i\end{align}$$

So I get a positive sign for $k=1$ and a negative sign for $k=3$.


Finding the edge pixels on only one side of the object

I found a way to detect the right hand side edge with MATLAB:

1.- acquiring image

A=imread('001.jpg');
figure(1);imshow(A)

A1=A(:,:,1);
[sz1 sz2]=size(A1)

figure(2);[P,h1]=imcontour(A1,1);

this detects the contour, but function imcontour produces an additional small amount of points completely unrelated to the contour, and quite far away

hold all;
plot(P(1,:),P(2,:),'r*')  % check

2.- quik-fix for imcontour producing a few points far away.

Px=P(1,:);Py=P(2,:);
nx_out=find(Px>sz2)
P(:,nx_out)=[]
ny_out=find(Py>sz1)
P(:,ny_out)=[]
hold all;
figure(3);plot(P(1,:),P(2,:),'r*');axis ij  % check

3.- generating image with the complete contour only and 1 pixel thick contour only

Px=P(1,:);Py=P(2,:);
A2=zeros(sz1,sz2)
P=floor(P)
for k=1:1:numel(Px)
    A2(Py(k),Px(k))=255;
end
figure(4);imshow(A2)

4.- Selecting right hand side edge only

P2=[0;0]
for k=sz1:-1:1
    L1=A2(k,:)

    if sum(L1)>0
        nL1=find(L1>0)
        P2=[P2 [k;nL1(end)]];
    end
end
P2(:,1)=[]

figure(5);
imshow(A2);
hold all;

for k=1:1:length(P2)
    plot(P2(2,k),P2(1,k),'r*')  % check
end

The result can be improved refining the way pixels are selected, instead of

P2=[P2 [k;nL1(end)]];

using

P2=[P2 [k;nL1(2:end)]];

but this is assuming that the excluded left edge is a fairly straight up line, that may not be the case.

For more complicated shapes, to first find the single pixel thick contour
 have a look at my answer: https://uk.mathworks.com/matlabcentral/answers/350136-how-t-oget-straight-edge-in-binary-image?s_tid=srchtitle


Help in interpreting the auto correlation graph

In my opinion this question is too broad. The first point is that you want to see if a given time series is random or not. I have doubt on how we can define this term here. Also you don't mention any distribution. Normally, we are interested to see if a time-series belongs to a specific distribution or not, for which there are standard techniques.

Regarding dependency, you should note that correlation only captures linear dependencies. You can have random variables that are highly dependent but have very low correlation. So it is important that what specific type of dependency you are looking for.

Edit: That said, referring to the ACF a peak at lag zero is visible. The reason is that at zero lag the time-series shows maximum similarity which is not surprising since it is obviously similar to itself! Note that the amplitude of this peak for a general complex signal $x[n]$ is $$R_{XX}[0]=\sum{x[n]x^*[n]}=\sum{|x[n]|^2}$$ which is the signal's power. You can easily show (theoretically) that $$R_{XX}[k]\le R_{XX}[0], \,\forall k$$

In addition, you should note that the autocorrelation of a white process (i.e. a process with uncorrelated random variables) is a pure delta function. So closer your autocorrelation to a delta means it is more uncorrelated, which is a good thing for a PRNG.

Regarding, your question about mutual information and entropy, yes both of them are used. Entropy is a well known measure for randomness. I cannot give you more details in here but can suggest you looking at this or similar resources.


Consider an IQ signal with amplitude or quadrature mismatch such as these QAM waveforms depicted below:

These can be decomposed into the sum of a perfectly balanced IQ signal and another signal that is modulating 0 to 180° with any fixed phase rotation. The perfect balanced IQ signal will be translated from the positive RF spectrum to baseband (as I depicted in the other post here Can Carrier Offset cause Image Problems) as a complex signal and thus if there was a carrier offset that would be maintained as an offset in the complex baseband signal but easy to correct for with a subsequent complex frequency translation. The signal that is on one axis only is a (potentially phase rotated) real signal and this will have alias distortion.

Consider an example with:

$$A\cos(\omega t)- jB\sin(\omega t))$$

If $A=B$ then this is simply $Ae^{-j\omega t}$. But if $B = (1+\Delta)A$ for example representing an amplitude imbalance of $\Delta$ then we get:

$$A\cos(\omega t)- j(1+\Delta)A\sin(\omega t))$$

$$=Ae^{-j\omega t} - j\Delta\sin(\omega t)$$

So the down-conversion that would have been image free with perfect IQ balance for a real input as follows:

$$r(t) = 2\cos(\omega_c t)$$
$$LO = \hat r(t) = e^{-j\omega_{LO} t}$$

$$r(t)\hat r(t) = 2\cos(\omega_c t)e^{-j\omega_{LO} t} = (e^{j\omega_c t}+e^{-j\omega_c t})e^{-j\omega_{LO} t}= e^{(\omega_c-\omega_{LO})t}+ e^{(-\omega_c-\omega_{LO})t}$$

Where the first component $e^{(\omega_c-\omega_{LO})t}$ is the only one that would appear at baseband.

With the amplitude impbalance as depicted above this example now becomes:

$$(e^{j\omega_c t}+e^{-j\omega_c t})(e^{-j\omega_{LO} t}- j\Delta\sin(\omega t)) $$

Consisting of the same result above with a new term given as:

$$(e^{j\omega_c t}+e^{-j\omega_c t})(-j\Delta\sin(\omega t))) = -(e^{j\omega_c t}+e^{-j\omega_c t})\frac{\Delta}{2}(e^{j\omega_{LO} t} - e^{-j\omega_{LO} t})$$

$$= -\frac{\Delta}{2}[e^{j(\omega_c+\omega_{LO})t}+ e^{j(\omega_c-\omega_{LO})t}+e^{j(-\omega_c+\omega_{LO})t}+e^{j(-\omega_c-\omega_{LO})t}]$$

Where we see we get two components at baseband, one $-\frac{\Delta}{2}e^{j(\omega_c-\omega_{LO})t}$ is constructive with our original signal (in this case decreasing the amplitude but no other distortion, while the other is an image $-\frac{\Delta}{2}e^{j(-\omega_c+\omega_{LO})t}$ which would cause serious distortion if the RF spectrums were not symmetric (as is typical for quadrature modulated signals) or if we had any frequency error between the carrier $\omega_c$ and the LO $\omega_{LO}$.

You can create a similar example with quadrature error, as well as the appearance of LO leakage with DC offset.

Consider this with the spectrums in a Zero-IF receiver as depicted below; the LO is depicted as a negative frequency only - with quadrature or amplitude error the positive frequency would start to appear starting with a very low level. You can follow the related spectrums when that occurs knowing the output is the convolution of the RF input spectrum with the LO spectrum to see the effect on the output.

With the super-heterodyne those overlapping images above would be separated appearing at their IF frequencies (and balanced since an IQ down-conversion would not be necessary from RF). But even a super-heterodyne receiver will have this same issue ultimately with quadrature and amplitude balance from the down-conversion at the IF frequency, since it does the down-conversion in two stages: first from RF to IF where we would then effectively have a Zero-IF down-converter from the IF frequency to baseband (either digitally or in the analog). However traditionally it has been a lot easier to achieve better amplitude and phase balance in the down-converter at the lower IF frequencies, as well as tighter filtering hence the superior performance of the super-heterodyne.

Thankfully with modern signal processing this is easy to detect and correct for, which has increased the popularity and use of Zero-IF receivers (and Direct RF transmitters). Below shows a typically correction approach for amplitude and phase imbalance:


Can you detect circles only using the hough line transform?

No this won't work. Simply because we do not know the center of the circle in advance (no oracle telling us) and therefore we cannot find a voting space parameterization where $r$ is fixed and only $\theta$ changes. For each point both quantities will change and this will still generate a curve in the Hough space. Having said that, I think below there is something interesting which uses tangent lines to detect circles:

Although what we will be doing is not exactly a Hough line transform, your way of thinking is on the right track. There is a strong connection there but because the tangent is not well represented in the Hough-circle transform, it probably won't work. Yet, what you have described brought me to the following idea:

It is possible to create an accumulator in the cartesian space ($X$-$Y$), where each pixel casts multiple votes as a partial line segment, along its gradient direction (or opposite gradient direction). The center of the circle will lie at the intersection of those segments. The width of the line segment can be related to the radius of the circle we like to detect. I think this approach is tightly coupled to operating in the $\textit{dual space}$ instead of the original and has demonstrated advantages, especially for conic fitting (well, a circle is basically a conic). A dual space is nothing but the space composed of the tangents. Each tangent is a point, just like in Hough.

Out of interest and because I think this is not a bad idea, I implemented this operation in MATLAB. Of course, instead of explicit edge-detection, I would simply use a gradient magnitude thresholding and use the image gradient as a cue of the direction of $\mathbf{NA}$ vector.

Here is the original test image:

And here is the accumulator space I get when I search for circles with $r=30$ pixels:

With peaks detected by the houghpeaks function:

The same peaks visualized back on the image:

It is possible to see that thanks to the voting along the line, the approach is pretty robust to the noise in the radius, i.e. we do not need to specify the radius explicitly.

For interested reader, I also provide the MATLAB code I used to generate those images:

function [ x,y,V ] = hough_circle_dual( I, maxRadius, threshold, numMaxCircles )

[M, N] = size(I);
[Gx, Gy] = imgradientxy(I);
[Gmag, Gdir] = imgradient(Gx, Gy);

V = zeros(M, N);

for i=1:M
    for j=1:N
        if (Gmag(i,j)>threshold)
            g = [Gx(i, j) Gy(i, j)]./Gmag(i,j);
            p2 = maxRadius.*g + [j, i];
            [x, y] = bresenham(j,i,p2(1), p2(2));
            ind = (x>0 & y>0 & x<=N & y<=M);
            x = x(ind); y = y(ind);
            ind = sub2ind([M,N], y, x);
            V(ind) = V(ind)+1;
        end
    end
end
peaks = houghpeaks(V, numMaxCircles);
x = peaks(:,1);
y = peaks(:,2);
end

% generates a bresenham line (compactly written)
function [x, y]=bresenham(x1,y1,x2,y2)
x1=round(x1); x2=round(x2);
y1=round(y1); y2=round(y2);
dx=abs(x2-x1); dy=abs(y2-y1);
steep=abs(dy)>abs(dx);
if steep t=dx;dx=dy;dy=t; end

if dy==0, q=zeros(dx+1,1);
else, q=[0;diff(mod([floor(dx/2):-dy:-dy*dx+floor(dx/2)]',dx))>=0];
end

if steep
    if y1<=y2 y=[y1:y2]'; else y=[y1:-1:y2]'; end
    if x1<=x2 x=x1+cumsum(q);else x=x1-cumsum(q); end
else
    if x1<=x2 x=[x1:x2]'; else x=[x1:-1:x2]'; end
    if y1<=y2 y=y1+cumsum(q);else y=y1-cumsum(q); end
end
end

Here is how I would call this function in MATLAB 2016:

I=imread('coins.png');
[x, y, V] = hough_circle_dual(I, 30, 10, 10);
figure,imshow(V,[]);
hold on, plot(y,x,'r+');

The Bresenham line generation used here is taken from the FileExchange post.


Will signal get distorted if channel bandwidth is lesser than the signal bandwidth?

The key here is the title of the book you are reading: "digital communicatinon by Taub schilling " as you say

What the book is talking about is digital communications and not analog communications in which the fidelity of the transmitted signal is of great importance; we want the output signal to be as close a replica of the input signal as is possible. In digital communications, it does not matter diddly-squat if the transmitted waveform is distorted as long as from the received waveform one can determine which of the two (more generally, few) possible transmitted waveforms $s_0(t)$ and $s_1(t)$ caused this received signal to occur. If the receiver can make this determination with high accuracy, then we have managed to transmit one bit (or more generally $\log_2$(few bits)) with high accuracy by choosing $s_0(t)$ or $s_1(t)$ according as we wish to send a $0$ or a $1$. Even if the received signal
$r(t)$ is badly distorted from what was sent, the receiver can still determine whether the transmitter sent $s_0(t)$ or $s_1(t)$: fidelity of transmission is unimportant.

With the above in mind, note that the Taub and Schilling description is dreadfully ambiguous, and in many spots downright incorrect. An equalizer (a concept from digital communications) cannot perfectly undo the distortion caused by forcing a 1 kHz analog signal through a low-pass filter with 3 dB bandwidth $1$ Hz. An inverse (analog) filter is trivial to write down --- it is a filter whose transfer function is $\frac{1}{H(f)}$ where $H(f)$ is the channel response --- but in most cases that inverse filter is unrealizable. The truth of the matter is that any filtering by a channel causes the transmitted signal to be spread out in time and reversing this time spreading and bringing the signal back to its pristine shape is in almost all cases impossible: once the genie is out of the lamp, he cannot be stuffed back in, ditto Pandora's box for those preferring a more Western analogy. That last sentence "The signal is then recoverable precisely as transmitted" is sheer nonsense.


An even function is symmetric:

$$f(x) = f(-x)$$

An odd function is anti-symmetric:

$$-f(x) = f(-x)$$

A signal in quadrature with another has all its Fourier series components phase-shifted by 90 degrees. The two signals are thus orthogonal.

To get the quadrature signal we apply the Hilbert transform (or Riesz transform in more than 1 dimension)

$$g(x) = \pm \mathcal{H}[f](x)$$

A cosine wave is even and symmetric. A sine wave is odd and anti-symmetric. Therefore if you apply the Hilbert transform to an even signal (cosine Fourier series components) you get an odd signal (sine Fourier series components). If you apply the Hilbert transform again you get the negative original signal. Apply again you get the negative of the original Hilbert transformed signal. Apply again you get the original signal. This is where the "quad" part comes in.

This is useful in signal analysis. Since the Hilbert transform only phase-shifts the Fourier series components, the energy of the signal remains constant. We can thus reconstruct the original signal (e.g. wavelet denoising). Also, given an even filter that responds to peaks and troughs, we can create an odd filter that responds to edges. The pair of filters is called a quadrature filter pair and allows the analysis of signal features. See the analytic signal for more details.


The answer is Non-uniform discrete Fourier transform

I suggest you to take a look in here: NFFT library. Tutorial for that purpose: NFFT 3.0 Tutorial.

You can also find a Python wrapper: pyNFFT.

  The generalisations of the NFFT include

  NNFFT - nonequispaced in time and frequency fast Fourier transform,

  NFCT/NFST - nonequispaced fast (co)sine transform,

  NSFFT - nonequispaced sparse fast Fourier transform,

  FPT - fast polynomial transform,

  NFSFT - nonequispaced fast spherical Fourier transform,

  NFSOFT - nonequispaced fast Fourier transform on the rotation group

  The NFFT is a C subroutine library for computing the nonequispaced
  discrete Fourier transform (NDFT) in one or more dimensions, of
  arbitrary input size, and of complex data. New: A Matlab interface is
  part of the NFFT3. We believe that our library, which is free
  software, and based on  FFTW (FFTW 3.x) should become the NFFT library
  of choice for most applications.


Your low-pass filter's pass band is smaller than your signal's bandwidth, so it is destroying a significant portion of your signal.  Given that, you cannot reliably reproduce the original signal without using sophisticated techniques like error correction codes.

I would try increasing $\alpha$ to 2400 Hz.  Then the rolloff point will be 2400 * .5 = 1200.  The resulting signal should look much better.  If that filter bandwidth is too large you can find your sweet spot somewhere between 1200 and 2400 Hz.


Comparison/Difference between two distinct plots

One way to compare similarly-shaped data that is translated in time is to look at the correlation between them.

If I do this in scilab:

ddata = diff(data',1,'c');

XC = xcorr(ddata(1,:),ddata(2,:));
CC1 = xcorr(ddata(1,:),ddata(1,:));
CC2 = xcorr(ddata(2,:),ddata(2,:));

den = sqrt(max(CC1)*max(CC2))

clf
plot(XC/den);
plot(CC1/max(CC1),'r')
plot(CC2/max(CC2),'m')

disp(max(XC)/den)

then I get this:

where the blue plot is the (normalized) cross-correlation of the time-differenced data and the red and magenta plots are the (normalized) auto-correlation of each individual time series.

The result is that the "correlation" between them is 0.8416782.


Confusion in initial condition of differential equation using Laplace transform transform

Initial conditions are always given at $t=0^-$, because they define the state of the system before any input is applied, and - by definition - the input is applied at $t=0$. The state at $t=0^+$ is determined by the initial conditions as well as by the input signal.
The unilateral Laplace transform can be used to solve LCCDEs with initial conditions $y(0^-), y'(0^-),\ldots$ because of the definition
$$\mathcal{L}\{f(t)\}=F(s)=\int_{0^{\color{red} -}}^{\infty}f(t)e^{-st}dt\tag{1}$$
from which it follows that
$$\mathcal{L}\{f'(t)\}=sF(s)-y(0^-)\tag{2}$$
Note that it's common to write initial conditions as $y(0),y'(0),\ldots$, when actually $t=0^{-}$ is meant.
EDIT: Concerning the example in the book: if there is no Dirac delta impulse in the current $i(t)$ at $t=0$, the capacitor voltage $v_c(t)$ cannot jump at $t=0$. Consequently, $v_c(0^-)=v_c(0)=v_c(0^+)$ must hold.
In general, the initial conditions define values of the output signal and its derivatives right before the source signal is switched on. It can be the case that the limits of the output and its derivatives exist at $t=0$. If that is the case, it doesn't make a difference if we use $t=0^{-}$ or $t=0$ or $t=0^{+}$ because the corresponding function values are all the same. If the output signal or its derivatives are discontinuous at $t=0$, then the distinction becomes important, and the initial conditions define the values just before the discontinuity, i.e., at $t=0^{-}$.


Deriving expression for the DTFT of a rectangular window

There's a mistake in the denominator of B. This is the way it should be done:

$$\begin{align}R(\omega)&=\frac{1-e^{-j\omega (M+1)}}{1-e^{-j\omega}}\\&=\frac{e^{-j\omega(M+1)/2}}{e^{-j\omega/2}}\cdot\frac{e^{j\omega(M+1)/2}-e^{-j\omega(M+1)/2}}{e^{j\omega/2}-e^{-j\omega/2}}\\&=e^{-j\omega M/2}\frac{\sin\left(\frac{M+1}{2}\omega\right)}{\sin(\omega/2)}\end{align}$$


I don't know if you are familiar with statistical signal processing and therefore will write my answer assuming that you are not. Everything I explain here is much better presented in any book about statistics. I would recommend Kay's book about detection theory. I first summarize your question by reformulating the 2 points you made, first in comprehensive English and then from a signal processing point of view.

You have access to two signals:

The signal you are looking for. Be it music, speech, or whatever. The only important assumption about this signal is that you perfectly know it. We will name it the signal of interest and denote it $s$. Let's say that it consists in $N$ samples denoted $s[n]$, with $0 \le n < N$.
A recording that contains the signal of interest as well as many other sources. The recording will thereafter be named the observation signal and denoted $x$. The other sources are all considered as noise in that they consist in anything but the signal of interest, and we will denote the noise $w$. The observation is longer than the signal of interest. Let's say it consists in $K$ samples denoted $x[k]$, with $K \ge N$ and $0 \le k < K$. Same goes for $w$.

You said you know what you are looking for ($s$), but you don't know  how badly it is "drown" ($w$) in what you are observing ($x$). From this point, the first thing I would try is a matched filter (MF) approach. I can't say for sure it will work for your problem because you gave us nothing about $w$. What I can say for sure is that this is the easiest approach that is very likely to work, both from a theoretical and a programming point of view. If this does not work, then you can come back here again. The fun part will start from there :)

But for now, let's focus on the matched filter. $s$ is shorter that $x$, so you should be able to take $M$ successive snapshots of $x$ of the same length as $s$, each denoted $x_m$. From $s$ and those $M$ snapshots, you can construct a third signal $y$ of length $M$, with each sample defined as follows :
$$y[m] = \sum_{n=0}^{N-1}{x_m[n]s[n]}$$

If $s$ is really in $x$, then you may see a clear maximum in $y$, at the index $m$ of the snapshot that contains $s$. If not, two possible conclusions:

$s$ is not in $x$. But that should not be the case, according to the presentation you gave of your problem.
Either the "level" of $s$ is way too low as compared to $w$ (we talk about a bad signal-to-noise ratio), or the statistical properties of $w$ are making this more complicated. The fun starts in the latter case.


Draw a frequency response graph from a simple equation

The tool of choice here is the z-transform, which is applied to your equation. It transforms sequences of complex numbers to functions in the complex plane. Its most useful properties are that it takes delays to powers of the argument z, and for z = exp(i omega) we get the fourier transform of the transformed sequence. Look up the details, I'll just show you how it works here.

y[n] = x[n] - x[n-486]  ---> Y(z) = X(z) - X(z)*z^(-486)

The equation in z-domain factors on the right hand side to give Y(z) = X(z) ( 1 - z^(-486) ).

Since Y(z) is the output and X(z) is the input, for linear systems we define the transfer function to be Y(z)/X(z) = H(z) and the above equation becomes

H(z) = Y(z)/X(z) = 1 - z^(-486)

Now if we substitute z = exp(i omega) we get the frequency response of your system, which is

H_f(omega) = 1 - exp(i omega)^(-486)

and with the properties of exp() it simplifies to

H_f(omega) = 1 - exp(- 486 i omega)

If you need actual frequency f instead angular frequency omega, just substitute omega = 2 pi f/fs, with sampling rate fs.

Finally, you'll probably want to know just the magnitude of the response, so look at abs(H_f).
In this specific case it's very simple to calculate because abs(H_f) = sqrt( H_f * conj(H_f) ) and conj(H_f) = 1 - exp(+ 486 i omega).

As requested, some code. Here's what you can do in matlab/octave to get the magnitude frequency response with a logarithmic frequency axis and magnitude axis in dB.

% Create a logarithmically spaced frequency axis with range 20Hz to 20kHz

fAxis = logspace( log10( 20 ) , log10( 20000) , 10000);

% calculate the complex frequency response at these frequencies

% Set SampleRate before calling this, like SampleRate = 44100;

Hf = 1 - exp(1i*486*2*pi*fAxis/SampleRate);

% plot the magnitude spectrum

semilogx(fAxis,10*log10(Hf.*conj(Hf)));


It can mean the following:
Consider any signal $x_{N}[n]$ of length $N$, and form the length $2N$ signal $x_{2N}[n]$ by repeating a copy of $x_{N}[n]$, then compute $2N$-Point DFT of this new signal $x_{2N}[n]$ whose every odd bin will be zero. i.e. $$X_{2N}[k] = 0 ~~,~~ \text{for}~~ k=1,3,5...,2N-1$$

whose demonstration can be performed with the following Octave/Matlab code:

 N = 16;
 x = randn(1,N);         % form an arbitrary signal x[n]
 x2 = [x x];             % repeat it once
 X2 = fft(x2,2*N);       % take DFT of length 2N
 figure,stem(abs(X2));   % display the result

Therefore it can mean that if you are performing an $M$ point DFT on your original signal $x[n]$, the signal $x[n]$ could actually be periodic with $\frac M2$ for $M$ even.

A more mathematical derivation of the result can be seen as follows:

Consider a signal $x_N[n]$ of length $N$ and its repetition $x_{2N}[n]$ of length $2N$ such that $x_{2N}[n]=x_N[n \mod N]$ for $n=0,1,\ldots,2N-1$.

Let's denote $N$-Point DFT of $x_N[n]$ as $X_N[k]$ for $k=0,1,\ldots,N-1$ and $2N$-Point DFT of $x_{2N}[n]$ as $X_{2N}[k]$ for $k=0,1,\ldots,2N-1$:

\begin{align}
&{\text{Start with defining $2N$-Point DFT of $x_{2N}[n]\ldots$}}\\
X_{2N}[k] &= \sum_{n=0}^{2N-1} {x_{2N}[n] e^{-j\frac{2\pi}{2N}nk}} \\ \\
&{\text{first split the sum and then substitude $n=n-N$ in the $2^{nd}$}}\\
&= \sum_{n=0}^{N-1} {x_{2N}[n] e^{-j\frac{2\pi}{N}n(k/2)}} + \sum_{n=N}^{2N-1} {x_{2N}[n] e^{-j\frac{2\pi}{2N}nk}} \\ \\
&{\text{recognise $x_{2N}[n]=x_N[n]$, $x_{2N}[n+N]=x_N[n]$ }}\\
&= \sum_{n=0}^{N-1} {x_{2N}[n] e^{-j\frac{2\pi}{N}n(k/2)}} + \sum_{n=0}^{N-1} {x_{2N}[n+N] e^{-j\frac{2\pi}{2N}(n+N)k}} \\ \\
&{\text{Expand the 2nd sum's multiplier factor}}\\
&= \sum_{n=0}^{N-1} {x_N[n] e^{-j\frac{2\pi}{N}n(k/2)}} + \sum_{n=0}^{N-1} {x_N[n] e^{-j\frac{2\pi}{N}(n+N)(k/2)}} \\ \\
&{\text{Replace $e^{-2j\frac\pi N}$ with $W_N$ for simplicity }} \\
&= \sum_{n=0}^{N-1} {x_N[n] W_N^{nk/2}} + W_N^{Nk/2}\sum_{n=0}^{N-1} {x_N[n] W_N^{nk/2}} \\ \\
&{\text{Recognise the sums as $N$-Point DFT $X[k]$ of $x[n]$, at $\frac k2$}}\\
&= X_N[k/2] + e^{-j\pi k} X_N[k/2] \\ \\
&= X_N[k/2] \cdot \big( 1 + (-1)^k \big)  \\ \\
&= \begin{cases}
2 X_N[k/2] , &\scriptstyle{\text{k=0,2,4,...,2N-2 }}\\
0  , &\scriptstyle{\text{k=1,3,5,...,2N-1 }}\\
\end{cases}
\end{align}

Shows that for such a defined signal , its odd DFT bins will be zero...


There are two issues here. First of all, you go from continuous time to discrete time, and second, you go from a continuous spectrum to a sampled spectrum. Note that the impulse response of an ideal discrete-time Hilbert transformer with frequency response

$$H(e^{j\omega})=-j\text{ sgn}(\omega),\qquad \omega\in[-\pi,\pi]\tag{1}$$

is given by

$$h[n]=\frac{1}{2\pi}\int_{-\pi}^{\pi}H(e^{j\omega})e^{jn\omega}d\omega=\begin{cases}\frac{2}{n\pi},& n\textrm{ odd}\\0,&n\textrm{ even}\end{cases}\tag{2}$$

So the zero samples at even indices don't come as a surprise.

However, with your method you don't get $(2)$ because you sample the spectrum, which results in aliasing in the time domain. You also get a non-zero imaginary part because your discrete spectrum does not satisfy the condition

$$X[k]=X^*[N-k]\tag{3}$$

where $^*$ denotes complex conjugation. Using the following definition for $X[k]$ will give you a purely real-valued time-domain sequence (up to numerical errors):

$$X[k]=\begin{cases}0,&n\in \{0,N/2\}\\
-j,& n=1,\ldots,N/2-1\\
j,&n=N/2+1,\ldots,N-1\end{cases}\tag{4}$$

where I've assumed that $N$ is even.

The inverse discrete Fourier transform (IDFT) of $(4)$ is given by

$$x[n]=\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j2\pi kn/N}=\begin{cases}\frac{2}{N}\cot\left(\frac{\pi n}{N}\right),&n\text{ odd}\\0,& n\text{ even}\end{cases}\tag{5}$$

Consequently, it is expression $(5)$ which you should expect instead of a sampled version of $1/\pi t$. Of course, $(5)$ approximates $(2)$, but the approximation gets worse for increasing magnitude of the index $n$ ($n\in [-N/2,N/2-1])$:


How to cover all spectrum in wavelet when it which is behalved?

It is covered - its just the signal that is the residual part of what you didn't filter yet.


Even though this may not be homework, it is a basic textbook problem, and for this reason I just give you a few hints, which should nevertheless help you solve the problem if you put in a bit of effort.

All questions refer to the filter's frequency response, which is simply obtained by evaluating $G(s)$ for $s=j\omega$.
The first question should read "Show that this property ... $\omega=0$, $\omega={\infty}$" (not $\omega=1$).
Your expression for $|G(j\omega)|$ is wrong; you wrote down $G(j\omega_0)$ (without magnitude). Of course, if $G(j\omega)=0$, or $G(j\omega)=1$, then also $|G(j\omega)|=0$ and $|G(j\omega)|=1$, respectively.
In order to see what happens for $\omega\rightarrow\infty$, write down $|G(j\omega)|^2$ (that's maybe easier than the magnitude), and take the limit.
For the 3 dB cut-off frequencies (there must be two positive solutions), use $|G(j\omega)|^2$ derived above and solve $|G(j\omega_c)|^2=\frac12$. Substitute $x=\omega^2$ and solve the resulting quadratic equation for $x$. From that solution, derive the two positive solutions for $\omega_c$.

All of this involves basic math skills, such as computing the (squared) magnitude of a complex number, and solving quadratic equations. If you find that difficult, read up on these basics because you'll need them all of the time when doing signal processing.


You can plot samples of a continuous density function just as you would do with any continuous function. Assuming $w(x)$ represents the continuous-time unit step function,

$$w(x-\Delta) = \begin{cases}{ 1 ~~~,~~~x \geq \Delta \\ 0 ~~~,~~~  \text{otherwise} } \end{cases}$$

then the shifted exponential pdf can be plotted with the following matlab/octave code :

clc; clear all; close all

K1 = -1;        % evaluation interval begins
K2 = 9;         % evaluation interval ends
N = 1000;       % number of points to display
x = linspace(K1,K2, N); % x = domain to use for displaying

lam = 1;        % exponential PDF parameter
delta = 2;      % shift amount
f = lam*exp(-lam*(x-delta));
f = f.*(x>delta);   % to implement w(x-delta)

figure, plot(x,f);
title('shifted exponential pdf for \lambda = 1, \Delta = 2');

The output will be:


KF is actually a mixture of a deterministic state propagator and a statistical estimator.

Despite it's name including the term filter, Kalman filter is not a simple frequency selective one. It's indeed a statistical recursive estimator of a state of a (linear) dynamic system. Yet on a broader sense it's called as a filter as it will separate a desired signal (the state) from corrupted measurements; hence filtering.

A linear dynamical system (LDS) is a physical construct whose behaviour is described by (linear) differential equations;

$$ \dot{X} = A X + Bu + W $$

where the vector $X$ is the state of the system, the matrix $A$ is the state propagation information (derived based on the differential equation of the LDS) , $u$ is a vector of deterministic inputs and $W$ is the state noise; depicting the uncertainity in the LDS propagation model $A$.

KF takes advantage of this deterministic behaviour of a LDS before estimating the current state $X$ purely from noisy measurements.

$$ Y = C X + Du + V $$

where $Y$ is the observation about the system and $V$ is the measurement noise on it.

Given a previous estimate of the state, it first deterministically propagates this previous estimate into current estimate through the matrix $A$ (process equation), and then performs a new finer estimate on the difference of the current measurement and propagated estimate using the measurment equation.


An easy way to deduce the nonzero range of discrete-time convolutions is the following:
Consider convolution of $x[n]$ with $y[n]$: $$z[n] = x[n] \star y[n] .$$
Let the shorther sequence $x[n]$ has a nonzero range of $M_1 \leq n \leq M_2$ , while the nonzero range of longer $y[n]$ be $N_1 \leq n \leq N_2$ .
1- Expand $x[n]$ as weighted impulses: $x[n] = \sum_{k=M_1}^{M_2} x[k] \delta[n-k]$ .
2- Distribute the convolution over expanded $x[n]$:
$$ (\delta[n-M_1]+...+\delta[n-M_2]) \star y[n] = \delta[n-M_1] \star y[n] + ... + \delta[n-M_2] \star y[n] $$
3- Consider the first impulse $\delta[n-M_1]$, and the last impulse $\delta[n-M_2]$.
4- Apply shifting property for those two impulses :
$$ \delta[n-M_1] \star y[n] = y[n-M_1]  ~~~,~~~  \delta[n-M_1] \star y[n] = y[n-M_2] .$$
5- Determine the nonzero ranges of $y[n-M_1]$ and $y[n-M_2]$:
$y[n]$ is nonzero over $N_1 \leq n \leq N_2$, then $y[n-M_1]$ will be nonzero for $$ N_1 \leq n-M_1 \leq N_2  ~~~~~~\Rightarrow ~~~~~~  N_1 + M_1 \leq n \leq N_2 + M_1 .$$
Similarly for $y[n-M_2]$ the nonzero range will be: $$ N_1 \leq n-M_2 \leq N_2  ~~~~~~\Rightarrow ~~~~~~ N_1 + M_2 \leq n \leq N_2 + M_2 .$$
One can see that the output $z[n]$ has a nonzero range that begins at the leftmost sample of $y[n-M_1]$, and ends at the rightmost sample of $y[n-M_2]$. Yielding:
$$N_1+M_1 \leq n \leq N_2 + M_2 .$$
Length of $z[n]$ can be found as: $$L_z = L_x + L_y -1 = N_2+M_2 - N_1-M_1 + 1 .$$
Applying this to yourcase, $N_1 = 0$, $M_1 = 0$, $N_2 = 3$ and $M_2=2$ which yields the nonzero range to be
$$ 0 + 0 = 0 \leq n \leq 5 = 2 + 3 .$$


Total Harmonic Distortion calculation and its origins

You can create non-linear digital systems (an example would be a system that finds the absolute value of the input). You can also simulate an analog non-linear system using DSP. The easiest way is to use a power series. Assuming that the non-linear system you're modeling is time-invariant and memoryless, you can approximate it as follows:

y(t) = a0 + a1*x(t) + a2*x(t).*x(t) + a3*x(t).*x(t).*x(t) + ....

where y(t) is the system output, x(t) is the input, and an are coefficients that depend on the system you're modeling.

As an example, you can do a simple simulation of mild harmonic distortion by setting x(t) equal to the a cosine signal, and setting a1=2, a3=-0.1, and all other coefficients equal to 0. a1 is the system gain, and a3 is the non-linear gain. By setting a3 to larger negative values, you increase the nonlinear effect -- you can think of it as an amplifier getting closer to saturation.

If you plot the spectra of x(t) and y(t), you'll clearly see harmonic distortion occuring.


there are issues.  given this convention for the continuous Fourier transform (and inverse)

$$ \mathscr{F} \Big\{ x(t) \Big\} \triangleq X(f) \triangleq \int\limits_{-\infty}^{+\infty} x(t) e^{-j 2 \pi f t} \ \mathrm{d}t $$

$$ \mathscr{F}^{-1} \Big\{ X(f) \Big\} \triangleq x(t) = \int\limits_{-\infty}^{+\infty} X(f) e^{+j 2 \pi f t} \ \mathrm{d}f $$

it changes the path of integration from the real axis to something else.  this comes up when using this fact:

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

to get, along with using scaling, this result:

$$ \mathscr{F} \Big\{ e^{j \pi t^2} \Big\} = \sqrt{j} \, e^{-j \pi f^2} $$

which is a linearly-swept "chirp" and it's spectrum.


How Exactly Does MATLAB Zero Pad Signal?

Lets say you have a vector $ x = {\left[ 1, 2, 3, 4 \right]}^{T} $.
You want to have a look on its DFT transform then you apply DFT on it and have the 4 points DFT transform of the data.
In MATLAB it would be something like:
vXDft = fft(vX);

Yet, what happens if you want the 64 point DFT of the same data?
In MATLAB it would be something like:
vXDft = fft(vX, 64);

So what happens?
MATLAB just adds zero elements at the end of the vector until its size is 64:
vXDft = fft([vX; zeros([60, 1])]);

This will match what happens above (I assumed Column Vector).
In practice you get a denser grid in the Frequency Domain (64 Points).
Where the points are interpolated by Dirichlet Kernel (The Discrete equivalent of the Sinc function).


Frequency Response Error for IIR Butterworth in GNU Octave

Don't know about GNU but in Matlab the frequency argument is specified relative to the Nyquist frequency (fs/2) and not the sample rate itself (fs). This makes the maxium legal argument for a relative cutoff frequency 1 (instead of 0.5). Your frequency argument should be [flo fhi]/fs*2.


What are the most robust implementations of passive/self-localization assuming high noise and reverberation?

It is well known that the GCC-PHAT algorithm (based on signals correlation) is robust in the case of reverberation. You can also use preprocessing algorithms to reduce noise.


Audacity Spectogram view: how to see exact frequency and intensity under the cursor?

What you are asking is possible using heuristics but impossible in general because of the uncertainty principle. A single sample is just a value, it does not even make sense to ask about its frequency or frequencies. Roughly speaking (though not strictly true), the concept of "frequency" necessarily implies some form of repetition, thus it only applies to a range of samples, not a single sample. You can, of course, look at a fixed-length window around the sample you are interested in, but the result (both the number of frequencies and their accuracy) will depend a lot on the window length.

I think Sonic Visualiser has pretty decent heuristics for instantaneous frequencies.


Cross-correlation or cross-covariance of non-zero mean signals

What are reasons to choose for cross-correlation or cross-covariance when comparing signals with non-zero mean?

Well, part of the issue is that cross-correlation as defined in your equation:

$$(f \star g)[n]\ \stackrel{\mathrm{def}}{=} \sum_{m=-\infty}^{\infty} f^*[m]\ g[m+n].$$

will not exist (or be infinite) if $f$ and $g$ have non-zero mean. So, in order to do the calculation, you really should subtract the mean beforehand.

If you're dealing with finite length signals, it's even worse, because you end up with an underlying triangular offset as per the image below.

This shows the sample cross-correlation between two gaussian processes with non-zero mean.  The black line is the sample cross-correlation calculated without subtracting the mean.  The red line is the sample cross-correlation calculated by first subtracting the mean from both signals.

  Does cross-correlation have any significance in this case, because $f\star g$ has mean $\mu_f\mu_g$?

It doesn't generally, unless you're very interested in the mean value.

And, as I said, as per your definition, any value of $\mu_f$ and $\mu_g$ both different from zero will result in an infinite sample cross-correlation.

  How would one find out the cross-correlation of signals that are not WSS, e.g. the mean is a linear trend? What would happen if we'd ignore this linear trend and simply subtract the signal means, then calculate cross-correlation? What if we subtracted the linear trend instead of the mean? What are the consequences for such mild violations of WSS-ity?

When looking at linear trends, it's usual to detrend before looking at the statistics of the variation around the trend. Sometimes doing this allows the variations to be treated as WSS.

If you keep the trend, you'll end up with a sample cross-correlation that looks like a parabola (the integration of a linear trend).  In this case, simply removing the mean doesn't help. See similar plot to first below, this time with two linear trend noisy signals.

Please see this exposition on the Mathworks site for a Matlab version of the above / below.

R Code Only Below

#34778
#install.packages('SynchWave')
require(SynchWave)

with_mean_ccf <- function(f,g)
{
  length <- length(f) + length(g) - 1

  ff <- c(f, rep(0,length - length(f) + 1))
  gg <- c(f, rep(0,length - length(g) + 1))

  FF <- fft(ff)
  GG <- fft(gg)

  ccf <- fft(FF * Conj(GG), inverse = TRUE)

  return(fftshift(Re(ccf)))
}

T <- 1000
mu_f <- 1
mu_g <- 2

f <- rnorm(T,1) + mu_f
g <- rnorm(T,1) + mu_g

with_mean <- with_mean_ccf(f,g)
no_mean <- with_mean_ccf(f-mean(f), g-mean(g))
plot(with_mean)
points(no_mean, col='red')
title('Sample crosscorrelation with (red) and without (black) mean subtraction')

a <- rnorm(T,1) + seq(1,T)/3
b <- rnorm(T,1) + seq(1,T)/2

with_mean_and_trend <- with_mean_ccf(a,b)
no_mean_and_trend <- with_mean_ccf(a-mean(a),b-mean(b))

plot(with_mean_and_trend, ylim=c(min(with_mean_and_trend,no_mean_and_trend), max(with_mean_and_trend,no_mean_and_trend)) )
points(no_mean_and_trend, col="red")
title('Sample cross-correlation of two linear trend signals')


Do you recognize this discrete realization of a low pass filter?

This is a discrete time filter and you have a block diagram, so forget the s-domain, and just work directly with difference equations and the Z transform.
In your block diagram, let's call the signal going into the delay block $w[n]$ and the signal coming out of the division block $v[n]$.
Writing some initial difference equations:
$$\begin{align*}
w[n] &= w[n-1] + v[n] \quad \Rightarrow \quad v[n] = w[n] - w[n-1]\\
\\
v[n] &= \dfrac{1}{k_1} u[n] - \dfrac{1}{k_1} w[n-1] \\
\\
y[n] &= w[n] - k_2 v[n]\\
\end{align*}$$
Performing some simple algebra, we get:
$$\begin{align*}
y[n] &= \left(1 - k_2\right) w[n] +k_2 w[n-1]\\
\\
w[n] - \left(1 - \dfrac{1}{k_1}\right)w[n-1] &= \dfrac{1}{k_1}u[n]\\
\end{align*}$$
So take the Z transforms and manipulate:
$$\begin{align*}
Y(z) &= \left(1 - k_2\right) W(z) +k_2 z^{-1}W(z)\\
\\
\dfrac{Y(z)}{W(z)} &= \left(1 - k_2\right)  +k_2 z^{-1} \\
\\
\dfrac{Y(z)}{W(z)} &= \left(1 - k_2\right) \left(1 - \dfrac{k_2}{k_2-1}z^{-1}\right)\\
\\
W(z) - \left(1-\dfrac{1}{k_1}\right)z^{-1} W(z) &= \dfrac{1}{k_1}U(z)\\
\\
\dfrac{W(z)}{U(z)} &= \dfrac{1}{k_1} \cdot \dfrac{1}{1-\dfrac{k_1 -1}{k_1}z^{-1} }\\
\end{align*}$$
The final transfer function is then:
$$\begin{align*}
\dfrac{Y(z)}{U(z)} &= \dfrac{Y(z)}{W(z)}\cdot\dfrac{W(z)}{U(z)}\\
\\
&= \dfrac{1-k_2}{k_1} \cdot \dfrac{1 - \dfrac{k_2}{k_2-1}z^{-1}}{1-\dfrac{k_1 -1}{k_1}z^{-1} }\\
\\
&= \dfrac{0.5}{0.5 + \frac{\tau}{T_s}}\cdot \dfrac{1 - (-1)z^{-1}}{1-\dfrac{\left(0.5 + \frac{\tau}{T_s}\right) -1}{0.5 + \frac{\tau}{T_s}}z^{-1} }\\
\end{align*}$$
You have a zero at the Nyquist frequency, $-1$ on the unit circle, so this is definitely looking like a lowpass filter.
You have a pole on the real axis just inside the unit circle, near $1$, at
$$ z_p = \dfrac{\left(0.5 + \frac{\tau}{T_s}\right) -1}{0.5 + \frac{\tau}{T_s}} = 1 - \dfrac{1}{0.5 + \frac{\tau}{T_s}} $$
which should provide a strong response near DC.  This is what you would expect from a lowpass filter.
You can do some algebra and take an inverse Z-transform to get the final difference equation:
$$\begin{align*}
\left( 1-\dfrac{\left(0.5 + \frac{\tau}{T_s}\right) -1}{0.5 + \frac{\tau}{T_s}}z^{-1}\right)Y(z) = \dfrac{0.5}{0.5 + \frac{\tau}{T_s}}\cdot \left(1 - (-1)z^{-1} \right) U(z) \\
\\
y[n] - \dfrac{\left(0.5 + \frac{\tau}{T_s}\right) -1}{0.5 + \frac{\tau}{T_s}} y[n-1] = \dfrac{0.5}{0.5 + \frac{\tau}{T_s}} u[n] + \dfrac{0.5}{0.5 + \frac{\tau}{T_s}} u[n-1]\\
\\
y[n] = \left(1-\dfrac{1}{0.5 + \frac{\tau}{T_s}}\right) y[n-1] + \dfrac{0.5}{0.5 + \frac{\tau}{T_s}} u[n] + \dfrac{0.5}{0.5 + \frac{\tau}{T_s}} u[n-1] \\
\end{align*}$$
So the filter is a scaled average of the two most recent input values with decaying feedback of the output value.
This could have been implemented with a much more straightforward block topology than the wonky one you were given.
In fact this difference equation is a 2 sample Moving Average filter, followed by an Exponentially Weighted Moving Average (EWMA) filter:
$$\begin{align*}
x[n] &= 0.5u[n] + 0.5u[n-1] \quad \text{(Moving Average)}\\
\\
y[n] &= \left(1-\dfrac{1}{0.5 + \frac{\tau}{T_s}}\right) y[n-1] + \dfrac{1}{0.5 + \frac{\tau}{T_s}} x[n] \quad \text{(EWMA)} \\
\end{align*}$$
From this answer: https://dsp.stackexchange.com/a/40465/28112 the -3 dB cutoff frequency of your EWMA filter is:
$$f_{3dB} = \dfrac{1}{2\pi T_s} \cos^{-1}\left[\dfrac{\left(\dfrac{1}{0.5 + \frac{\tau}{T_s}}\right)^2}{2\left(1-\dfrac{1}{0.5 + \frac{\tau}{T_s}}\right)}\right]$$
Two other observations for you:
$k_2$ sets the balance of the weight given to the inputs $u[n]$ vs. $u[n-1]$
$k_1$ sets the $\alpha$ of the EWMA filter, which sets the balance of the weight given to the previous output vs. the new input.
Update
So now that we know where the digital filter's pole is, let's work through from the continuous time filter's pole, to see if it matches, using the bilinear transform relationship.
$$\begin{align*}
s_p = -\dfrac{1}{\tau} &= \dfrac{2}{T_s} \cdot \dfrac{z_p -1}{z_p+1}\\
\\
z_p+1 &= -2\dfrac{\tau}{T_s}\left(z_p-1\right) \\
\\
z_p\left(1 + 2\dfrac{\tau}{T_s}\right) &= -1 + 2\dfrac{\tau}{T_s} \\
\\
z_p\left(1 + 2\dfrac{\tau}{T_s}\right) &= -2 + \left(1 + 2\dfrac{\tau}{T_s}\right) \\
\\
z_p\left(0.5 + \dfrac{\tau}{T_s}\right) &= \left(0.5 + \dfrac{\tau}{T_s}\right) -1 \\
\\
z_p & = \dfrac{\left(0.5 + \dfrac{\tau}{T_s}\right) -1}{0.5 + \dfrac{\tau}{T_s}} = \dfrac{k_1 -1}{k_1}
\end{align*}$$
So that checks out. Someone did apparently start from a single pole filter of the form
$$H_a(s) = g\cdot\dfrac{1}{1+\tau s}$$
$H_a(s)$ has a zero as $s \rightarrow \infty$ (along any direction from the origin of the s-plane).  Looking again at the bilinear transform relationship:
$$\begin{align*}
s_0 &= \dfrac{2}{T_s} \cdot \dfrac{z_0 -1}{z_0+1}\\
\\
\lim_{z_0 \to -1} s_0 &\rightarrow -\infty\\
\\
z_0 &= -1\\
\end{align*}$$
Which also agrees with the zero we found analyzing the digital filter.
Update 2
And just to give some insight as to why $k_2 = 0.5$
$$\begin{align*}
s_0 &= \dfrac{2}{T_s} \cdot \dfrac{z_0 -1}{z_0 +1}\\
\\
&= \dfrac{2}{T_s} \cdot \dfrac{\frac{k_2}{k_2-1} -1}{\frac{k_2}{k_2-1} +1} \\
&= \dfrac{2}{T_s} \cdot \dfrac{k_2 - k_2 + 1}{k_2 + k_2-1}\\
\\
&= \dfrac{2}{T_s} \cdot \dfrac{ 1}{2k_2 - 1}\\
\\
&= \dfrac{1}{T_s} \cdot \dfrac{ 1}{k_2 - 0.5}\\
\end{align*}$$


log-likelihhood function for N sample of data

I’m assuming that the trace is what is confusing, so let me give you some clues.

First an identity $$ Trace(\mathbf{ABC)})= Trace(\mathbf{BCA}) $$

and for scaler $x$,
$$ Trace(x) =x$$

Let $\mathbf{a}$ and $\mathbf{b}$ be $N \times 1$ vectors and $\mathbf{C} $is a $N \times N$ matrix, so $$ \mathbf{ a^T C b} = Trace ( \mathbf
{C \,b a^T}) $$

The rest you should be able to figure out using the definition of a Gaussian, independence, and properties of the log.


OCR (optical character recognition) for Indian languages - segmentation

Why do not use hough transform for finding lines and then finding table region?

you can use hough transform to find horizontal and vertical lines. and then extract region of lines.


As mentioned in another answer, image erosion is quite efficient to fill the gaps between the wall hatched lines. The key point is to choose a structure element that is big enough to allow connection to be done.

You should use a structure element with both horizontal and vertical dimensions of at least half the distance between hatched lines (try different values!). You may end with a final erosion with a small structure element to remove remaining noise (the thin and lonely lines in your case).

Then it is your responsibility to detect and refine the borders ;)

Here is attached the result I obtained with a square structure element of 19x19 pixels on a sample of the image you linked (logical OR one) assuming that the distance between lines is 19 pixels:

Hint: The output is smoother with a larger structure element but artefacts appear at the corners.


One possibility is use something similar to an MFCC Mel triangular filterbank with your FFT input vector, but using log spacing instead of Mel spacing for the triangular filter centers.  See: http://cmusphinx.sourceforge.net/doc/sphinx4/edu/cmu/sphinx/frontend/frequencywarp/MelFrequencyFilterBank.html and https://www.mathworks.com/matlabcentral/answers/195651-creating-mel-triangular-filters-function, and modify the triangle frequencies as needed.


In the simplest terms, the slowness or fastness of a filter refers to its time domain characterization. Specifically a fast filter will have a short impulse response so that the output goes to steady state as quick as possible in time, whereas a slow filter will have a long impulse response, whose steady state output occurs after a longer delay, so that it will not be responsive to transient changes in the input.

Now according to the time-bandwidth product feature of LTI filters (which is based on a Fourier transform property) if a filter has wide frequency support, then its time support will be narrow, and vice versa. As a consequence, low cutoff frequency filters will have narrow frequency support and hence large time support. But a large time support is equivalent to long impulse response which is stated to be a slow filter as above.

Therefore a low cutoff filter (a narrow frequency support) will be a slow responding one, whereas a larger cutoff filter will be fast responding one.


Getting HRTF from 3-D scan + acoustic simulation

There are various school of thoughts about this. A German guy named Klaus Genuit came up with a parametric model of HRTFs that's based on the actual physcial geometry of head, shoulder, pinna, etc. He actually now runs a company that craeates binaural products and measurement gear, see for example
http://www.head-acoustics.de/downloads/publications/binaural_technology/DAGA_99_Physical_Modeling_of_individual_HRTFs.pdf

Based on this work you could potentially create individual HRTFs using a some 3D scanning scheme. I do vaguely remember a company coming doing exactly that for a bunch of people at my employer at the time. I didn't think the HRTFs that came out that way were particularly good but I didn't do a lot of work with them either.

Individualized HRTFs only help to a certain extent. Let's say you want to record a concert. You can do this binaurally using a dummy head. However, there is currently no way to convert this recording to your individual HRTFs after the fact, even if you have them.


Term for integral multiples of fundamental period?

If $T$ is a positive real number such that
$$x(t) = x(t+T) ~\text{for all}~t \in \mathbb R\tag{1}$$
for continuous-time signals, or if $T$ is a positive integer and
$$x[n] = x[n+T]~\text{for all}~n \in \mathbb Z \tag{2}$$ for discrete-time signals, then $x$ is said to be periodic with period $T$ or to have period $T$.
A signal $x$ that has period $T$ also has period $2T$, period $3T$, $\ldots$, but might (or might not) have smaller period such as $T/2$ or $T/3$ etc. The smallest value of $T$ for which $(1)$ or $(2)$ is called the fundamental period of the signal, and all other periods necessarily are integer multiples of the fundamental period.

So, to answer your question about the name for integer multiples of the fundamental period, they can be called a period of $x$ which is perfectly valid, or, if it is necessary to emphasize the fact that $T$ is not the fundamental period, one can say more prolixly, "$T$ is a period but not the fundamental period of $x$." The more one want to nitpick (or strive for clarity and completeness), the more wordy one must be.

Note that if $T$ is the fundamental period of a signal $x$ and if $T^\prime = mT$ (where $m$ is an integer greater than $1$) is a period, then we can express the continuous-time signal $x(t)$ of period  $T^\prime$ as a Fourier series with fundamental frequency $$f_0^\prime = \frac{1}{T^\prime} = \frac{1}{mT}= \frac{1}{m}f_0$$ where $f_0 = \frac 1T$ is the fundamental frequency of $x(t)$ of fundamental period $T$. $f_0^\prime$ is called a subharmonic of $f_0$.  Notice, however that the Fourier series that we get by regarding $x(t)$ as a signal of period $T^\prime$ is exactly the same as the Fourier series that we get by regarding $x(t)$ as a signal of period $T$. In
$$x(t) = \sum_{n=-\infty}^\infty c_n \exp(j2\pi nf_0^\prime t)
= \sum_{n=-\infty}^\infty c_n \exp(j2\pi \frac nm f_0t)\tag{3}$$
the $c_n$'s have value $0$ except when $n$ is an integer multiple $km$ of $m$ and so the right side of $(3)$ can be expressed as
$$x(t) = \sum_{k=-\infty}^\infty c_{km} \exp(j2\pi k f_0t)$$
which is precisely the Fourier series for $x(t)$ regarded as a signal of period $T$. In short, while we have introduced the concept of subharmonics of $f_0$, $x(t)$ has no frequency content at these subharmonic frequencies.

A similar result holds for discrete-time signals $x[n]$ but I am too lazy to write it all out. Perhaps @OlliNiemitalo will resuscitate his answer and spell it all out for those the benefit of those who are interested.


Firstly, because the speed can be positive and negative (meaning approaching to and going from sensor if target is moving), you need to shift your doppler indices from [0, N-1] to
[-N/2,N/2 - 1] by simple subtraction by N. Then you just multiply scale [-N/2,N/2-1] by:

lambda/(2*N*T_rpi)

where:
- lambda is wavelength for your FMCW radar (lambda = c/Fc, c - speed of light and Fc central transmit frequency of chirp - for 77-81GHz you would use 79GHz)
- N is number of chirps transmitted
- T_rpi is chirp repetion interval (meaning T_active + T_idle if you have idle time after chirp transmission)

If you use multiple transmitters in time division multiplexed mode, you would need to multiple by:

lambda/(2*N*T_rpi*n_tx)

You then gain scale in m/s.


Sending complex signal over a single wire

Note that signals sent over wires (and over the air, and over any medium) are always real. What quadrature means is that, on a passband channel (wired or otherwise), you can transmit two signals at the same time. The first signal, which we'll call $s_I(t)$, is mixed with a carrier $c_I(t)=\cos(2 \pi f_c t + \phi)$; the second signal, $s_Q(t)$, is mixed with $c_Q(t)=\cos(2 \pi f_c t + \phi + \pi/2)$, where $f_c$ is the center frequency of the channel.

Note that it is important to use a passband channel. Quadrature communication is impossible on the baseband channel. Maybe this is what confused you, since one tends to associate wires with baseband. However, as long as you don't exceed the bandwidth of the wire, it's perfectly possible to transmit modulated signals over them. In fact, in telephony, the very first trunk links used single sideband modulation to multiplex several phone calls over a single wire -- this was before PCM and digital telephony were even invented.

The transmitter would then calculate the following signal: $$s(t) = s_I(t)c_I(t) - s_Q(t)c_Q(t),$$ where the negative sign is just a convention. Assuming the carrier phase $\phi=0$ for simplicity, this signal can also be written as $$s(t)=\text{Re}[(s_I(t)+js_Q(t))e^{j 2 \pi f_c t})].$$ This expression is sometimes useful in a digital implementation, since you only have to keep track of two complex signals instead of four real signals.

The channel adds noise, so the received signal is $r(t) = s(t) + n(t)$. The noise is usually modeled as Gaussian, white noise with PSD equal to $N_0/2$.

In the receiver, we have (ommitting a lot of algebra readily found in textbooks and on other answers in this site): $$\hat{s}_I(t) = s_I(t) + n_I(t) = \text{LPF}[r(t) c_I(t)]$$ and $$\hat{s}_Q(t) = s_Q(t) + n_Q(t) = \text{LPF}[r(t) c_Q(t)],$$ which means that the receiver can recover both transmitted signals (albeit with noise).

How to actually get all this to work in an actual system requires solving other problems that have not been addressed here and that can take a whole textbook to explain. I recommend that you start with this free textbook: http://sethares.engr.wisc.edu/telebreak.html (link at the end of the page).


Hilbert-Huang starts with empirical mode decomposition (EMD). I know one HHT code is available on Matlab central but I personally find it not very robust and extremely sensitive to edge effects. A more robust solution may be found here : http://perso.ens-lyon.fr/patrick.flandrin/emd.html

Once the EMD is performed, the missing step to complete HHT is the Hilbert transform which is fairly straight forward.

More details and links are mentioned in this post.


What is the meaning of the length of diagonal lines in FFT

If you rotate an image (2D function), its 2D Fourier transform also rotates.

Those lines mean you have something like a barcode oriented in that direction in your image. In one direction it has lots of variations and behaves like a white noise with flat spectrum and in the other direction it has very slow variations.

I'm not sure what do you mean by the length of diagonals, but it's a visual measure of maximum frequency available in those directions.


invfreqz() and frequency domain filter specification

The warning is due to modeling a unit delay which has wide (constant) frequency response thus will have very large number differences (ill conditioned matrices) in the translation between domains and looking for an IIR solution to the FIR problem (note that the algorithm converges to a trivial but valid solution where the denominator coefficients match the numerator-- if the OP set $m$ to be larger, the algorithm will converge to populate those coefficients as well in the numerator). The recommended approach is to start with FIR only, and only introduce IIR if needed. Regardless the result does match the expected very well in both time and frequency as plotted below after running the OP's code:

Frequency Response

Impulse Response

The above is the IIR equivalent to the OP's FIR target given $m=2$ in the call to invfreqz, which would then converge to an IIR solution. To force an FIR solution, set $m=0$ which results in the following for this case:

b = 0 1 0 0 0 0 0
a = 1

(Note the zeros for b were actually +/-3e-17 or lower)

This is a useful approach for fitting to a precise frequency response curve. For general lowpass, bandpass and multiband filter design, I prefer to use the least squares firls() function with simple targets of 1 and 0 for the passbands and stopbands. However if a more elaborate response is desired (in a curve fitting sense), this is the approach I would take:

Start with a very large $n$ and $m=0$ and run invfreqz() to evaluate the resulting impulse response (b coefficients). Assuming the impulse response converges sufficiently to what can be assumed 0 (finite) at the start and end of the filter (evaluate in dB scale as small values can have a significant impact), reduce $n$ to this length and rerun invfrewqz. Evaluate the frequency response using freqz and increase $n$ if needed to meet target (or reduce $n$ until desired maximum error is reached). The result of this is the best FIR only solution that can be achieved. Given the coefficients of the FIR filter is the impulse response, any filter with an impulse response that converges to zero can be realized as an FIR only filter, but implementing as an IIR can reduce the total number of coefficients needed (with all the cost of IIR implementation such as possible instability, round-off error accumulation, dead-beat responses, etc).

If the number of elements in the FIR is unacceptably long such that an IIR solution is desired, increase $m$ and repeat which will result in a smaller $n$ and an IIR solution. (Make $m$ excessively long to start and evaluate b and a on a log scale to determine the number of significant coefficients). If using an IIR solution, it is recommended to then factor the resulting solution into 2nd order sections to minimize round-off errors that can accumulate significantly in higher order IIR filters. (Resulting in 2nd order biquad filters).


Remove noise introduced while converting a Delta-Sigma Modulated (DSM) signal into a pulse-code modulated (PCM) signal?

The spectrum appears to be the noise-shaping associated with Delta-Sigma Modulation. Specifically we see that the noise increase follows a 20 dB/decade slope, indicative of a first-order delta-sigma modulator.
The high end of the noise rolls off consistent with the high frequency roll-off of the filter applied.
Below are some charts I have explaining the first-order Delta-Sigma modulator and how such noise shaping occurs, here specifically with a Delta-Sigma Digital to Analog Converter (but it could easily be a 1 bit digital output just the same).

And equivalently the Delta Sigma Analog to Digital converter diagram is shown below, which will be used to develop a simple loop model for providing an understanding of how noise shaping occurs.

Below is a qualitative loop model, providing an intuitive understanding of noise shaping if one has an appreciation of filtering in simple control loop models. Specifically it is showing an 1 bit A/D Converter as a gain with quantization noise added (this is a highly non-linear process, therefore the model is qualitative versus a quantitative prediction).

The path through the loop from the quantization noise input ($N_q$) to the output follows a high pass function as shown below. Quantization noise would typically be white noise (uniform over all frequency) for a traditional data converter, while with the implementations using Delta Sigma modulation we get this resulting noise shaping resulting in significantly higher precision at a given sampling rate after the high frequency noise is filtered out.

Typically higher order implementations are used, such as the 2nd order Delta Sigma shown below (and higher orders such as 5th order are commonly implemented). This provides much steeper noise shaping (much less noise in the band of interest after filtering) and also serves to significantly reduce pattern noise which would manifest as spurs in the spectrum.


What you might be forgetting is that dwt does downsampling. After filtering (low-pass or high-pass), the filtered signal is subsampled. In order to keep the number of samples. So aliasing may occur, depending on the quality of the filters.

And db2 is a quite poor filter. So if you idwt either approx or details, by replacing the other coefficients by zeros, you almost get your coefficients in order: on the top, the original signal,  the two low frequencies in the second plot, the high frequencies in the third plot. Of course you have some aliasing, due to the poor resolution.

Try the attached code with db24, you will get much less artifact.

What happens now with the approximation and the details? My interpretation is that the approximation is now sampled at $500$ kHz, and the details too, except that it is now aliased to the base-band:

and they need to be unmixed by reciprocal filters of the idwt.

The FFTR.m code should be downloaded.

But why study sines with discrete  wavelets?

f = [20/1000  100/1000  300/1000];
nSample = 1024; time = (1:nSample)';
data = cos(2*pi*time*f);
dataSum = sum(data,2);
waveName = 'db2';
[a,d] = dwt(dataSum,waveName);
z = zeros(size(a));

dataSumD = idwt(z,d,waveName);
dataSumA = idwt(a,z,waveName);

timeAD = 1:2*length(a);
timeAD=timeAD(1:2:end);

figure(1)
subplot(3,2,1)
plot(time,dataSum);axis tight
subplot(3,2,2)
plot(FFTR(dataSum));axis tight

subplot(3,2,3)
plot(time,dataSumA);axis tight
subplot(3,2,4)
plot(FFTR(dataSumA));axis tight

subplot(3,2,5)
plot(time,dataSumD);axis tight
subplot(3,2,6)
plot(FFTR(dataSumD));axis tight

figure(2)
subplot(2,2,1)
plot(timeAD,a);axis tight
subplot(2,2,2)
plot(FFTR(a,2));axis tight

subplot(2,2,3)
plot(timeAD,d);axis tight
subplot(2,2,4)
plot(FFTR(d,2));axis tight

So you have the code, I do not know whether you did stuff correctly, and wavelets are quite useful if used widely on signals that are worth  it. Sines typically are not on the  list.


Can we get the same result from between partial FFT?

I don't think "separable" means what you think it means – it means that you can implement the 2 DFT as the concatenation of a $N_{rows}$ amount of $N_{columns}$ long DFTs followed by $N_{columns}$ DFTs of length $N_{rows}$.

So, no, your approach is not correct.


Behaviour of the poles in transfer function $H(s)$, given system properties

A ROC can be chosen to make the transfer function causal and/or stable depending on the pole/zero plot. If the ROC extends outward from the outermost pole, then the system is causal.


Transmitting over Unknown Fading Channel

The strategy depends heavily your "realistic scenario", e.g. which estimator, which equalizer, which estimation error, etc.

In general scenario, you cannot trust the estimate of $h$, it means that you cannot use coheent detection. Thus try non-coherent detection (detection by energy) : to send bit $0$, use 2 channel uses (2 symbols) $x_0 = [x[0], x[1]] = [1, 0]$ and to send bit $1$, use  $x_1 = [x[0], x[1]] = [0, 1]$. The receive symbols are $y = [y[0], y[1]]$

In the well-known Rayleigh flat fading channel, i.e. $y[m] = h[m]x[m]+w[m]$, the optimal detection rule will be
\begin{align}
\hat{x} &= x_0 \textrm{ if } |y[0]|^2 > |y[1]|^2\\
\hat{x} &= x_1 \textrm{ otherwise}
\end{align}


Why DTFT coefficients are periodic and why continuous Fourier transform coefficients are not periodic?

The explanation by @Maximilian Matthé is a standard and formal approach for this question. But I think it is not intuitive and easy to understand the reason. In the following, I will try to explain from another aspect.

First of all, periodicity means infinity. For a signal in the temporal/spatial domain, the periodicity is somewhat easy for realization, such as extending it periodically and infinitely. But for the frequency domain, it is quite different. Different position means different component with different changing speed. Therefore, if periodicity is expected in the frequency domain, infinitely rapid changing is expected in the temporal/spatial domain. In the real world of our life, most signals changes smoothly, thus periodicity is impossible in the frequency domain. While for discrete-time signals, the model is based on the $\delta(t)$ function, for which the changing is infinitely sharp and rapid, thus periodicity becomes possible in this case.


Maximum Likelihood Detection of Signal Vectors in Gaussian Noise

The probability of error of the ML detector is equal to the probability that the received vector $\mathbf{y}$ is closer to $\mathbf{x}_1$ than to $\mathbf{x}_0$, which is equal to the probability that the noise component in the direction of $\mathbf{x}_0-\mathbf{x}_1$ is greater than half of the Euclidean distance between $\mathbf{x}_0$ and $\mathbf{x}_1$. The Euclidean distance between $\mathbf{x}_0$ and $\mathbf{x}_1$ is

$$D=||\mathbf{x}_0-\mathbf{x}_1||=\sqrt{\sum_{i=1}^d(x_{0,i}-x_{1,i})^2}=\sqrt{\sum_{i=1}^d4E_s}=2\sqrt{dE_s}\tag{1}$$

The noise variance in the direction of $\mathbf{x}_0-\mathbf{x}_1$ equals $N_0/2$ (as in any other direction). The probability that a zero mean Gaussian noise variable with variance $N_0/2$ assumes a value greater than $D/2$ is given by

$$P_E=Q\left(\frac{D/2}{\sqrt{N_0/2}}\right)=Q\left(\frac{\sqrt{dE_s}}{\sqrt{N_0/2}}\right)=Q\left(\sqrt{\frac{2dE_s}{N_0}}\right)\tag{2}$$


Basis signals for the discrete-time Fourier Series

The complex exponentials

$$\phi_k[n]=e^{j\frac{2\pi k}{N}n},\qquad k=0,1,\ldots, N-1\tag{1}$$

form an orthogonal basis, i.e.

$$\sum_{n=0}^{N-1}\phi_k[n]\phi_l^*[n]=0,\qquad k\neq l\tag{2}$$

So clearly they are all different, and they are even orthogonal (and hence not linearly dependent).

It's straightforward to show that Eq. $(2)$ is true:

$$\sum_{n=0}^{N-1}\phi_k[n]\phi_l^*[n]=\sum_{n=0}^{N-1}e^{j\frac{2\pi (k-l)}{N}n}=
\begin{cases}
\sum_{n=0}^{N-1}e^{j\frac{2\pi (0)}{N}n}=\sum_{n=0}^{N-1}1=N,\qquad & k = l \\
\\
\frac{1-e^{j2\pi (k-l)}}{1-e^{j\frac{2\pi (k-l)}{N}}}=\frac{1-1}{1-e^{j\frac{2\pi (k-l)}{N}}}=0,\qquad & k\neq l \\
\end{cases}\tag{3}$$


It is easy to see that using the trigonometric identity in Equation $(1)$ below
$$
\cos(\theta_1 \pm \theta_2) = \cos(\theta_1)\cos(\theta_2) \mp \sin(\theta_1)\sin(\theta_2) \tag{1}
$$
We have
\begin{align}
f(t)\cos\big(\pi(t - 1)\big) &= f(t)\overbrace{\cos(\pi t - \pi)}^{\text{use Equation}\ (1)}\\
&\equiv -f(t)\cos(\pi t)\tag{2}
\end{align}
Which can be rewritten as
$$
-f(t)\cos(\pi t) = -f(t)\cos\left(2\pi \frac t2\right) = -f(t)\cos(2\pi f_c t)\quad\text{where}\quad f_c = \frac 12
$$
And we know the Fourier cosine modulation frequency-shift property in equation $(3)$:
$$
\mathcal F\big\{x(t)\cos(2\pi f_c t)\big\} = \frac 12\big[X\left(f - f_c\right) + X\left(f + f_c\right)\big]\tag{3}
$$
You can now easily use Equation $(2)$ and $(3)$ to find the solution.


Discrete time Final Value Theorem applied to feedback system

I would derive the total transfer function directly in the transform domain. Your input-output equation can be written as
$$Y(z)\big(1-z^{-1}\big)=\alpha G(z)z^{-1}Y(z)+\beta z^{-1}X(z)\tag{1}$$
where $G(z)$ is the transfer function of the high-pass filter. From $(1)$ you directly obtain the transfer function
$$H(z)=\frac{Y(z)}{X(z)}=\frac{\beta z^{-1}}{1-z^{-1}\big(1+\alpha G(z)\big)}\tag{2}$$
which clearly has a pole at $z=1$ because $G(1)=0$.

The mistake in your calculation is in step 2. You increased the lower bound of the summation from $1$ to $2$ without changing the power of $\phi$ from $i$ to $i-1$. In the transfer function this results in a $\phi^2$ in the last part of the denominator instead of a $\phi$. The denominator should be:
$$D(z)=1 - (1+\Psi\phi)z^{-1}-\Psi(\phi-1)\frac{\phi z^{-2}}{1-\phi z^{-1}}\tag{3}$$
For $z=1$ this evaluates to
$$D(1)=1 - (1+\Psi\phi)-\Psi(\phi-1)\frac{\phi }{1-\phi }=1 - (1+\Psi\phi)+\Psi\phi=0\tag{4}$$


Program to display amplitude of sine signal at each interval and display separately in MATLAB

it's actually easy, you only need to find the maximum value of your sine, which is in turn happens to be your amplitude, then you can compare them with if an if statement.

this gives you the amplitude

max(signal);

use if statement to compare the two values

if max(signal1) == max(signal2)
%% your condition comes here
end

Hope this helps.

cheers.


To find coordinate or index of drew region, first you have to find the coordinate of 4 colored dots, which is easily obtained by using the transformPointsForward function over the intrinsic coordinate of these 4 points in your original image. intrinsic coordinate is the index of your image matrix.

then you have to do some geometric calculations. first, find the length of one of red lines which are the (length of full warped image)/(sin(theta) + cos(theta)), then distance of red line to the nearest colored point is (length of red line)*cos(theta)*sin(theta). theta is the rotation angle. using this distance and the coordinate of 4 colored points you could find the coordinate of your ROI corners.


By separately changing both the leading edge and the trailing edge of each pulse, you have converted the square wave into what can be viewed as an asymmetric double-edge pulse-width-modulated (DEPWM) signal. Analysis of the Fourier spectra of various types of deterministic PWM signals (including the type mentioned in the previous sentence) can be found in the paper Z. Song and D. V. Sarwate, "The Frequency Spectrum of Pulse Width Modulated Signals", Signal Processing, vol. 83, pp. 2227-2258, 2003. (The article is available behind the publisher's paywall but there are some free copies floating around the web).  The results show that the spectra of PWM signals do have the spikes shown in the OP's diagram just as the square wave spectra do, but there is additional clutter between the peaks. That is, the PWM Fourier spectrum is not a Fourier series like the square wave spectrum is; PWM spectra are of a continuous-time signal that contains periodic components (giving those spikes) as well as non-periodic components giving the clutter.  The details (and the nuances) are too messy to transcribe here; interested readers should read the original.


How to characterize this whitening filter?

Hi: The material you are asking about is on page 7 of the document at the link below.

http://www.ee.ic.ac.uk/hp/staff/dmb/courses/DSPDF/DSPDF_2010.pdf


Yes, you will need to make a 50Hz notch filter.

I would suggest you to look at this post.

If you tell what type of program you use (Python, R, etc) I could try to provide you some code aswell.


Standard deviation is varying with change in amplitude of DC signal

I believe your intuition is correct.

The filter you're using is a low pass filter with a DC gain of 1. Without seeing the output of the filtering operation, I suspect that you're seeing the effect of initial conditions in the filter: getting from 0 to 1 in your first instance takes fewer samples than getting from 0 to 10 in your second.

To test this, try taking the middle $L/2$ of your filtered signal and see what std(y(L/4:3*L/4) is then.


Isolating Bat Species without prior training base

1) Based on the above, does would this be an ideal solution in order
  to isolate the bats and therefore monitor their activity or would I
  need to implement further algorithms in order to isolate the bats
  calls?

The short answer is no, this would be no where near ideal. This type of work strongly overlaps with speech/voice recognition, and the algorithms (or more so, the features) involved are not trivial.

There are certainly frequency features that one can extract from audio signals of species to try and cluster them, and the DFT would certainly be the transform of choice for this type of raw feature. However the DFT will only give you frequency scores, although it where one would start. Other frequency metrics and features can be more complex. (Look at cepstral analysis, MFCC, etc).

  2) I have only ever had experience of using the Cooley-Tukey FFT
  algorithm, would this be sufficient enough for this type of project?

Remember that the DFT is a Transform, whilst the FFT is an efficient algorithm for this transform. Furthermore, the Cooley-Tukey Algorithm is one way of implementing an FFT. In your work, you will use canned libraries that have the FFT. You don't need to implement an FFT in the 21st century, just as you would not be expected to implement the square root function.

I strongly suspect that you will need to compute the Short-Time-Fourier-Transform (STFT) for this type of project. In anticipation of your next needs, please see the wiki link on the STFT, as well as this post here. The STFT allows you to analyse signals in the dual time-frequency domain. You will most likely want to start your feature extraction in that space.

  let's assume I have computed the STFT and the MFCC and extracted the
  features, how would I therefore identify the species of the bats?
  Would I therefore need some kind of training in order to complete this
  task?

Once you have decided on what sets of features you need to use, you will of course need a learning algorithm to discriminate them properly. Since you do not have a training set, there are two options. The first obvious is to attain a training set, and/or meticulously label your data based on what you know to be true. Once you have your labels, then you can use supervised training techniques like logistic regression or a simple perceptron to come up with separating hyperplanes in feature space that disjoint your data. The second option however, would be for you to use unsupervised learning techniques (such as K-means for example).

If you use something like K-means, you will need to know apriori how many different bat species you are trying to discriminate against. This is the value of your hyper-parameter $k$. The output of K-means will be $k$ vectors, each of them being in $N$-dimensional points in feature space, representing your cluster centroids. Now, every time you get a new audio sample, you:

Compute $N$ features and place this into a feature vector.
Compute the L2 vector norm, (aka, the Euclidean distance) between the given feature vector, and the $k$ centroid vectors you pre-computed from using $k$-means.
Pick the $i$th centroid out of the $k$ cetroids that corresponds to the closest distance to your point.

Mathematically, you want to declare a bat of type $i$, among $k$ type bats by:

$$
{arg\,min}_{i \in (1, 2, ...k) }  \sqrt{\sum_{n = 1}^{N} (\mathbf{x} - \mathbf{f_i})^2 }
$$

where $\mathbf{x}$ is an $N$ x $1$ feature vector of your unclassified audio sample, and $\mathbf{f_i}$ is the $i$th $N$ x $1$ cluster centroid found by $k$-means.


Determining the mean and standard deviation in real time

There's a flaw in Jason R's answer, which is discussed in Knuth's "Art of Computer Programming" vol. 2. The problem comes if you have a standard deviation which is a small fraction of the mean: the calculation of E(x^2) - (E(x)^2) suffers from severe sensitivity to floating point rounding errors.

You can even try this yourself in a Python script:

ofs = 1e9
A = [ofs+x for x in [1,-1,2,3,0,4.02,5]]
A2 = [x*x for x in A]
(sum(A2)/len(A))-(sum(A)/len(A))**2

I get -128.0 as an answer, which clearly isn't computationally valid, since the math predicts that the result should be nonnegative.

Knuth cites an approach (I don't remember the name of the inventor) for calculating running mean and standard deviation which goes something like this:

 initialize:
    m = 0;
    S = 0;
    n = 0;

 for each incoming sample x:
    prev_mean = m;
    n = n + 1;
    m = m + (x-m)/n;
    S = S + (x-m)*(x-prev_mean);

and then after each step, the value of m is the mean, and the standard deviation can be calculated as sqrt(S/n) or sqrt(S/n-1) depending on which is your favorite definition of standard deviation.

The equation I write above is slightly different than the one in Knuth, but it's computationally equivalent.

When I have a few more minutes, I'll code up the above formula in Python and show that you'll get a nonnegative answer (that hopefully is close to the correct value).

update: here it is.

test1.py:

import math

def stats(x):
  n = 0
  S = 0.0
  m = 0.0
  for x_i in x:
    n = n + 1
    m_prev = m
    m = m + (x_i - m) / n
    S = S + (x_i - m) * (x_i - m_prev)
  return {'mean': m, 'variance': S/n}

def naive_stats(x):
  S1 = sum(x)
  n = len(x)
  S2 = sum([x_i**2 for x_i in x])
  return {'mean': S1/n, 'variance': (S2/n - (S1/n)**2) }

x1 = [1,-1,2,3,0,4.02,5]
x2 = [x+1e9 for x in x1]

print "naive_stats:"
print naive_stats(x1)
print naive_stats(x2)

print "stats:"
print stats(x1)
print stats(x2)

result:

naive_stats:
{'variance': 4.0114775510204073, 'mean': 2.0028571428571427}
{'variance': -128.0, 'mean': 1000000002.0028572}
stats:
{'variance': 4.0114775510204073, 'mean': 2.0028571428571431}
{'variance': 4.0114775868357446, 'mean': 1000000002.0028571}

You'll note that there's still some rounding error, but it's not bad, whereas naive_stats just pukes.

edit: Just noticed Belisarius's comment citing Wikipedia which does mention the Knuth algorithm.


Why is DTFT of $e^{jn\omega_0}$ an impulse train?

First of all, note that the discrete-time Fourier transform (DTFT) of a sequence is always periodic with period $2\pi$:

$$X(\omega)=\sum_{n=-\infty}^{\infty}x[n]e^{-jn\omega}\tag{1}$$

This is obvious because $e^{-jn\omega}$ is $2\pi$-periodic in $\omega$.

Next, consider the inverse DTFT:

$$x[n]=\frac{1}{2\pi}\int_{-\pi}^{\pi}X(\omega)e^{jn\omega}d\omega\tag{2}$$

For $X(\omega)=\delta(\omega-\omega_0)$, we get from $(2)$

$$x[n]=\frac{1}{2\pi}\int_{-\pi}^{\pi}\delta(\omega-\omega_0)e^{jn\omega}d\omega=e^{jn\omega_0}\tag{3}$$

[Note that without loss of generality we assume that $-\pi\le\omega_0\le\pi$ because $X(\omega)$ is $2\pi$-periodic. $X(\omega)=\delta(\omega-\omega_0)$ is a shorthand for $X(\omega)=\sum_k\delta(\omega-\omega_0-2k\pi)$. That shorthand can be used here because we know that $X(\omega)$ is the DTFT of some sequence, and hence $2\pi$-periodicity is understood].

Since $\textrm{IDFT}\{\delta(\omega-\omega_0)\}=e^{jn\omega_0}$ we're inclined to accept

$$\textrm{DTFT}\{e^{jn\omega_0}\}=\sum_{n=-\infty}^{\infty}e^{jn\omega_0}e^{-jn\omega}=\delta(\omega-\omega_0)\tag{4}$$

(where $2\pi$-periodic continuation is understood), even though
the sum

$$\sum_{n=-\infty}^{\infty}e^{jn\omega_0}e^{-jn\omega}\tag{5}$$

doesn't converge in the conventional sense. However, this can be made mathematically sound by interpreting the sum $(5)$ as a distribution.


Why discretize a continuous transition matrix in Kalman Filter?

Have you taken a look at the documentation - section 2.2.4 discusses a Linear Kalman filter model that is very similar to the one you described. From that example you see that:

The resulting $\bf{A}$ matrix does not depend on the process noise
The $\bf{A}$ matrix only depends on the size of the time step.
The $\bf{Q}$ matrix only depends on the size of the time step, and the value of the continuous noise covariance, ie. diag(q) in the example

You can use discrete noise kinematic models, they are probably a bit easier to understand. The continuous noise model assume that the noise is applied throughout the time step interval i.e. from time $t_0$ to $t_1$. The discrete model assumes the noise is applied instantaneously at times $t_0$ and $t_1$ (and so on).

See chapter 6 -"Estimation for Kinematic Models" in the book - "Estimation with Applications to Tracking and Navigation" (Bar-Shalom, Rong Li, Kirubarajan) for a derivation of the continuous and discrete time versions of the state transition models.


Camera Calibration - why do you have to move the calibration board?

This is due to the optimization problem being rather high-dimensional (around 11 parameters). With only a single observation of the calibration board, there would be multiple possible combinations of parameters explaining the observed feature point locations (unless a very constrained camera-model is used). Only a sufficient number of sufficiently independent observations will properly constrain the optimization to yield a well-defined answer. Otherwise the routine will converge to some (local) minimum.

Think of this as solving a set of equations with more unknowns than equations. You will be able to find a solution, however, there would be an infinite number of solutions (one or more free parameters).

With respect to a board viewed from the front (frontoparallel orientation), filling the entire image. Now -- is this a calibration board close to the camera viewed with a wide-angle lens (short focal length), or is it a calibration board far away, viewed with a zoom lens (long focal length)? This shows that we need images with the board tilted (foreshortening).

With some constraints, you can calibrate a camera with a single image of two calibration patterns forming e.g. a 90 degree angle, or two perfectly parallel boards at differents depths.

I have written an article on the usual optimization process, which might help you understand how the process could end up not converging to the global solution: https://calib.io/blogs/knowledge-base/camera-calibration

This has little to do with the fixed grid of the imaging sensor as suggested in other answers.


The transpose (or more generally, the Hermitian or conjugate transpose) of a filter is simply the matched filter. This is found by time reversing the kernel and taking the conjugate of all the values.

Chapter 5 of Jon F. Claerbout - Earth Soundings Analysis: Processing versus Inversion has a great discussion of the topic of Adjoint operators (or transpose, or back projection, or matched filtering). It's in the context of Seismic imaging, but it really gets to the heart of what you're actually trying to do (in general) with these techniques.


How to develop a spectrogram (2D array) from audio signal?

Spectrogram is obtained using Short Time Fourier Transform technique. It is used to visualize frequency changes over time. Imagine you have a very long signal whose frequency content vary with time. You can be tempted to take the FFT of this whole data. But that will not show you the variation of frequency content over time. So you will have to split your signal into smaller chunks and take FFT of each chunk. That way you will be able to see how the frequency content varies over time. This is how the 2-D array is generated. For each chunk, you get the FFT output as one single column vector(in each column as per your figure). For consecutive chunks, you place them adjacent to each other. so as you move across row, you are moving across time.

How do you decide how to split your time signal into smaller chunks? Suppose your FFT size is $N$, you can split your whole signal into chunks of $N$ and take FFT of each without overlapping these chunks. So at time $n$, your FFT output is
$$
X(k,n) = \sum_{m=n}^{m=n+N}x[m]e^{-j2\pi km/N}
$$
When taking FFT, we have not done any pre-processing on $x[m]$. We just did a rectangular windowing of $x[m]$ from $m=n$ to $m=n+N$ before taking FFT. But rectangular windowing is not a good windowing technique because it will increase spectral leakages. So instead of that, you will apply other windowing techniques like Hann or Hamming or Kaiser to mitigate this spectral leakage. Windowing is just point-by-point multiplication of $x[m]$. For window function $w[n]$ , $0 \le n \le N-1$,
$$
X(k,n) = \sum_{m=n}^{m=n+N}x[m]w[m-n]e^{-j2\pi km/N}
$$
For the picture you showed,
$$
\text{Column, X[:,m]} \rightarrow  \text{ Spectral Content at time index m}\\
\text{Row, X[k,:]} \rightarrow \text{ Variation of frequency index k across time index m}\\
$$
Size of windowed $x[m]$ may be less than FFT size $N$ but for simplicity I am assuming same here.

There is one more parameter while obtaining this 2-D array. I mentioned earlier that the chunks do not overlap. But if the frequency content varies fast across each chunk, you may end up getting a jagged display when you plot the spectrogram. So it is imperative to provide some amount of overlap between chunks. The overlap size is $L$. So for each time index (row), the increment would be $L/T_s$, where $T_s$ is the sampling interval. To summarize, any two consecutive columns of the 2-D spectrogram array will be computed as
$$
X_r(k,n+rL) = \sum_{m=n+rL}^{m=n+rL+N}x[m]w[m-n-rL]e^{-j2\pi km/N}\\
X_{r+1}(k,n+(r+1)L)= \sum_{m=n+(r+1)L}^{m=n+(r+1)L+N}x[m]w[m-n-(r+1)L]e^{-j2\pi km/N}
$$


1) The first equation should look like: h_pred = ifft2 ( fft2(k) ./ fft2(x) ). You have a small typo there, I believe.

Make sure you first zero-pad the kernel to the size of image.

2) MATLAB also has a blind deconvolution function:
http://www.mathworks.com/help/images/ref/deconvblind.html

I don't know if you are referencing to this one, but for Toeplitz approach, I would refer to:

[1] P. C. Hansen, “Deconvolution and Regularization with Toeplitz Matrices,” Numerical Algorithms, vol. 29, no. 4, pp. 323–378, 2002.


Let's view this from a slightly different angle:

If the signal is properly sampled, i.e. in accordance with the Shannon/Nyquist criterion, then the samples contain all information about the original signal. If not, all bets are off, so we'll skip this for now.
Interpolation is then equivalent to sampling the signal at a non-integer time. In your case you want to know what is x(t=2.1) by using the information x[1], x[2], x[3] ...
There is only one "correct" answer to "what is x(t=2.1). In your example it happens to be 1.05061339.... It can be calculated using the Whittaker Shannon interpolation (see Wikipedia) as $$x(t)=\sum_{n=-\infty }^{\infty }x[n]\cdot \mbox{sinc}(\frac{t-nT}{T})$$
This can be interpreted as the convolution with an impulse response. However, it's not "finite".
Unfortunately The Whitaker Shannon interpolation is not very practical: The impulse response is infinite, the number of zeros is infinite and it's falloff with time is slow, so it's not a great candidate for windowing to get a time limited impulse response.
The Whittaker Shannon interpolation is equivalent to convolution with the impulse response of an ideal low pass filter. All practical interpolation methods will also involve  a low pass filter. Even linear interpolation can be interpreted as a low pass filter, it's just a very bad one.
Interpolation filters are made "practical" by adjusting filter parameters to the requirements of the application: how much precision do you need, what's the spectral content of the signal, what latency can you tolerate, etc.
Filter parameters to adjust are filter order, cutoff frequency, rolloff, filter design method (equiripple, least squares, etc), linear or minimum phase, side lobe height, stop band attenuation, pass band ripple, etc.


There are a few potential issues:

Your code doesn't compute a full output vector, you just keep the last three samples (often referred to as "state variables". Your output should be a vector of the same length as the input
This is only a second order filter. Depending on the exact properties of your signal, this  may not be enough to isolate the component that you want. I'd recommended a Butterworth bandpass filter. In this case you can independently control the lower cutoff and the higher cutoff of your bandpass and also the steepness of the slopes. There are butterworth design tools in "Octave", see http://octave.sourceforge.net/signal/function/butter.html
If you need to go to higher order filters, they should be implemented as "cascaded biquads" or "cascaded second order sections". There are tons of code examples to be found through google


Calculating the extent of spectral leakage

Windowing in the time domain results in a convolution with the transform of that window in the frequency domain.  So you need to compute the result of this convolution to see the effect (which some call "leakage").

You can sample the transform of the window function used (rectangular, Von Hann, etc.) to determine the FFT frequency response of a windowed pure single frequency sinusoid.  For a rectangular window, the magnitude of the apparent "leakage" in an FFT result will be samples of a normalized Sinc function (more precisely a pair of normalized periodic Sincs or Dirichlet kernels), sampled at integer offsets from the frequency bin center offset, multiplied by the sinusoid's magnitude.

An attempt at the equation for the full Dirichlet pair frequency response result is on my dsp web page here: http://www.nicholson.com/rhn/dsp.html#4

But using just a normalized Sinc as the window transform equation is usually close enough, except for the bins right next to DC or 0 Hz and to bin N/2.

For windows other that rectangular, use their transform instead of a Sinc function, up-sampled and/or interpolated if needed.


Unit impulse response of a cascade interconnection of three discrete-time systems

HINT:

Consider the cascade interconnection of three causal LTI systems, [...]"

What does that imply for the impulse response $h_1[n]$?
What is more, since the total impulse response is causal, and you know from its definition that $h_2[n]$ is causal, $h_1[n]$ must also be causal.


Why the RMS of a PSD curve is the root of the area below

I think this is simply an aspect of Parseval's Theorem (e.g. click me)

It simply says: sum of squares in the time time domain equals sum of squares in the frequency domain. Substitute "sum" for "integral" if using the analog domain. In other words: total energy in the time domain equals total energy in the frequency domain.

This can easily reproduce your formula, $g_{RMS}$ represents the  time domain energy and the integral on the right represents the energy in the frequency domain. The exact scaling depends on the details like length of the signal, periodicity, sample rate, etc.