How do I implement the Structure tensor (in Matlab)?

However the input to the Gaussian might have a mean different from 0.

Unless your image is completely constant, $I_x^2$ and $I_y^2$ are guaranteed to have a mean > 0 (assuming the input pixels are all reals, of course).

  How do I compensate for this?

Why would you want to do that?

  I was thinking substracting the mean first.

Subtracting the mean before applying a derivative operator like Sobel or Scharr? That shouldn't make any difference.

  The Scharr operators are supposed to be perfectly rotational symmetric

Perfect rotational symmetry? Is that even possible with a finite kernel? I think the Scharr operator is just optimal (in a least-squares sense) for that filter size. With larger filter sizes, you get better rotational symmetry.

  How do I Implement this easily in Matlab?

You can use imfilter to apply the Sobel or Scharr filters, then .* to multiply/square the results pixelwise, then imfilter again to apply the gaussian smoothing kernel.

  How about in C++?

Any decent image processing library should contain functions for linear filters and pixel-wise multiplication. e.g. for the IPP, the functions would be called ippiFilterRow/Column..., ippiMul... and ippiSqr...

  Any comments regarding choice of smoothing kernel and differentiation operator?

I think the gaussian smoothing kernel is pretty standard. Sobel and Scharr are both common choices for the differentiation operator. Another good choice would be a gaussian derivative filter. Which one's best depends on your application: how noisy are your images? How sharp are the contrasts you're looking for? (If there's lots of noise, and the contrasts are sharp, a derivative filter with a larger aperture like the gaussian derivative might be better). How important is rotational symmetry?

  How do I normalize the Gaussian derivative?

The gaussian is: $\frac{1}{2 \pi  \sigma ^2}e^{-\frac{x^2+y^2}{2 \sigma ^2}}$

Deriving by x yields:

$\partial _x\left(\frac{1}{2 \pi  \sigma ^2}e^{-\frac{x^2+y^2}{2 \sigma ^2}}\right) = -\frac{e^{-\frac{x^2+y^2}{2 \sigma ^2}} x}{2 \pi  \sigma ^4}$

So I think you should divide by gd = -x.*g/(si^2); because sum(g)=2*pi*(si^2)


Coordinates should be passed in y,x order, not x,y.
Full updated code:

#include "opencv2/core/core.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"

#include <iostream>

using namespace cv;
using namespace std;

Mat world;
Point startPt(0,0);
Point endPt(0,0);
float dist;

void onMouse( int event, int x, int y, int flags, void* )
{
    if( event == CV_EVENT_LBUTTONUP) startPt = Point(x,y);
    if( event == CV_EVENT_RBUTTONUP) {
        endPt   = Point(x,y);
        Vec3f s = world.at<Vec3f>(startPt.y, startPt.x);
        Vec3f e = world.at<Vec3f>(endPt.y, endPt.x);
        float dx = e[0]-s[0];
        float dy = e[1]-s[1];
        float dz = e[2]-s[2];
        dist = sqrt(dx*dx + dy*dy + dz*dz);
    }
}

int main( /*int argc, char* argv[]*/ ){
    VideoCapture capture;
    capture.open(CV_CAP_OPENNI);
    capture.set( CV_CAP_PROP_OPENNI_REGISTRATION , 0);
    unsigned t0=clock();
    if( !capture.isOpened() ){
        cout << "Can not open a capture object." << endl;
        return -1;
    }
    unsigned elapsed=clock()-t0;
    cout << "initialized in "<< elapsed <<" s. ready!" << endl;
    namedWindow( "depth", 1 );
    setMouseCallback( "depth", onMouse, 0 );
    for(;;){
        if( !capture.grab() ){
            cout << "Can not grab images." << endl;
            return -1;
        }else{
            Mat depthMap,show;
            capture.retrieve(world, CV_CAP_OPENNI_POINT_CLOUD_MAP);
            if( capture.retrieve( depthMap, CV_CAP_OPENNI_DEPTH_MAP ) ) depthMap.convertTo( show, CV_8UC1, 0.05f);
            line(show,startPt,endPt,Scalar(255));
            putText(show,format("distance: %f m",dist),Point(5,15),FONT_HERSHEY_PLAIN,1,Scalar(255));
            imshow("depth",show);
        }
        if( waitKey( 30 ) >= 0 )    break;
    }
}


How to distinguish between FFT power leak or close frequencies

Assuming the spectrum is sufficiently interpolated to see the detail within the resolution bandwidth (which is a reasonable assumption given the question), and the waveform isn't further windowed beyond the rectangular window (as I would argue that any such windowing that is applied is equivalent to the presence of multiple tones closely spaced given windowing is an applied amplitude modulation) then one obvious way to quickly distinguish these two cases is by observing the width of the main lobe (as well as all other features such as null locations, or non-existence of the nulls, but the main lobe width would be robust under noise conditions)-- given the kernel of the rectangular window (the DFT of the window ) convolves in frequency with each tone. If there is one tone we would see the kernel alone. If there are two tones spaced closer to together, then both tones would convolve with the kernel, and therefore widen the result.  (The kernel for a rectangular window in time is the Dirichlet Kernel in frequency, basically an aliased Sinc function, so if the spectral peaks are the same width as the Dirichlet Kernel, which is known, not computed, then we can easily declare that peak to be a single tone).
If sufficiently interpolated and noise free then we can also quickly detect this by the lack of frequency nulls (although the nulls are soon buried by any noise so this would not be a robust test). If the spectrum wasn't interpolated and the frequencies are so close that a bin width increase is not obvious, then we could compare the phase and magnitudes of the leakage values in all other bins to what would be expected with the underlying Dirichlet Kernel, under noise conditions we could use all bins and a correlation to what would be expected for a single tone to provide a confidence of one tone or not.
Further Details on using Main Lobe Width Measurement
This process of main lobe width measurement would be precise in confirming a single tone down to very small frequency offsets in the case of a single exponential tones of the form $e^{j\omega t}$, but will have a lower offset limit for the case of real tones due to the interaction of the two exponential tones that each real tone contains ($cos(\omega t) = 0.5(e^{j\omega t} + e^{-j\omega t}))$. In which case we would see even for a single tone a small variation in main lobe width versus frequency that becomes more pronounced as we approach the Nyquist boundaries. The effects of this, and interaction from any other tones further away for the same reason, become significantly reduced if the waveform was windowed prior to computing the DFT (which can't be assumed based on what the OP has stated, and as I already mentioned I would argue the result of windowing a single tone results in the presence of multiple closely spaced tones!). To show the practical extent of using the main lobe width approach, I plot the variability of the main lobe width versus frequency for the rectangular window case at the bottom of the post, which provides insight into the practical limitation of how close we can distinguish multiple tones (pretty close! down to 0.02 bins for most of the frequency range, and 0.2 bins for all frequencies except the end points).
To demonstrate this effect of simply measuring the width of the main lobe for the case of an interpolated spectrum, and its relative insensitivity to noise to address hotpaws comment, see this example below with a rectangular window and two equally leveled tones that are within one bin width (the resolution bandwidth of the rectangular window is 1 bin, so more than that would not be sufficient for "close"):

Here is the comparison of the composite signal with and without noise, where the noise is only 6 dB lower:

(And in this case, even if the spectrum was not interpolated, meaning the only samples available are at the integer bin locations; we can see how the increased width of the main lobe would still be obvious through inspection)
For further details on the limitations of the main-lobe width measurement for a real tone, the following plot below shows the variation of the main-lobe bin width at -3 dB for a rectangular windowed 128 point DFT (interpolated by zero-padding out to 12800 points), where the "truth" given by a single exponential tone approaches 0.886 bins for large $N$ (where $N$ is the number of samples, derived from $2\omega/\pi$ for $Sinc(\omega) = \sin(\omega)/\omega = 0.707$, and can be exactly determined for any size N from the Dirichlet Kernel, $D(\omega)$, from $2\omega/\pi$ for $D(\omega) = 0.707$, where:
$$D(\omega) = \frac{\sin(N\omega/2)}{N\sin(\omega/2)}$$
and $\omega$ is the normalized radian frequency in units of radians/sample : $\omega \in [0, 2\pi)$ for $N \in [0, N)$
With the results shown in the plot immediately below:

As given by the convolution, this bin width will be increased by the frequency separation of two tones that are less than 1 bin apart, and as long as that separation is greater than twice the variability of a single tone, it would be detectable (for example the variability when we are more than 10 bins away from the Nyquist boundaries is less than 0.01 bins, which holds as N is increased). This plot below was simulated but could be analytically predicted by the interference pattern of the underlying Sinc function that becomes the Dirichlet Kernel through aliasing. This demonstrates the practicality of using this for distinguishing a single tone from multiple tones closely spaced.


Zero-phase-shift filters must be symmetric about the origin

Consider the Inverse Fourier Transform as the impulse response of the filter. In the case of an ideal brick-wall filter with zero-phase, the impulse response would be a Sinc function centered about zero in time (given by the Inverse Fourier Transform).

This is non-causal (time domain response is greater than 0 for $t<0$), but in order to be zero-phase as shown, the time domain response must be complex conjugate symmetric.  It need not be real; in order for the frequency response to be real, the time response real components must be even and the imaginary components must be odd (complex conjugate symmetric).

But to visualize the above plot further and what "zero-phase" means, a delay in time is a linear phase in frequency, thus if we shift the impulse response to the right (as done prior to truncation for a causal solution), the magnitude response is unchanged but the phase in no longer "zero-phase". If we continued to see the next plot I have in this series, which shows the impulse response truncated by the dashed window in this plot below, we would then see that the magnitude response of the filter would also be distorted, but that is beyond the question being asked here.

I summarize the properties of odd and even functions as well as other universal Fourier Transform properties in the figure below. The time and frequency domains are interchangeable for these properties; for example, if a function is periodic in one domain, it is discrete in the other. If a function has only real components in one domain (zero-phase), then all real components in the other domain must be an even function about 0, and all imaginary components must be odd about zero (inverse sign). This results in complex conjugate symmetry.

The relationship of the Hilbert Transform mentioned for the causal and anticausal functions applies to minimum phase systems.

What is interesting is how these tie together and why I choose to describe the above relationships in terms of even and odd functions rather than simply saying the other domain will be complex conjugate symmetric. This way we can see how a causal function in time (or a single sided function in frequency) is the sum of even and odd functions, and therefore MUST be complex in the other domain given that we need both real and imaginary components to create such even and odd functions.

Regarding the OP’s second question: a notch filter usually implies a very narrow band of rejection to remove certain frequency components. For more details see this post here: Transfer function of second order notch filter


Your first version of $J(a*b)[n]$ is correct:

$$J(a*b)[n] = \sum_{p=0}^{N-1} a[p]b[-n-p]$$

because you simply have to flip $n$, and not the summation index.

The problem is in your formulation of the circular convolution of $(Ja)[n]$ and $(Jb)[n]$. It should be

$$((Ja)*(Jb))[n]=\sum_{p=0}^{N-1} a[-p]b[-n+p]$$

because you have to negate the arguments of both sequences, so $p$ becomes $-p$, and $n-p$ becomes $-n+p$. Now you can simply change the sign of $p$ because it doesn't matter in which order you compute the sum. This gives

$$((Ja)*(Jb))[n]=\sum_{p=0}^{N-1} a[p]b[-n-p]$$

which is identical to $J(a*b)[n]$ as given above, as it should be.


Producing masks from mostly-static-background subtraction when given the background

You can try working off a different representation of the image, other than color. You can diff the image gradient for example. It is very unlikely that the gradient is the same on the face as the background. A gradient based descriptor such as HOG or SIFT captured at a suitable scale over the image has worked well for me in the past. You can even combine different representations of the image, such as RGB value as well as HOG. IF the image is different in either mark as a change. You will still get the odd few pixels that don't come up as a change but you can smooth over the result with a smoothing kernel with a similar scale tot he size of the object you want to detect.


Are the 4 types of linear phase FIR filters the only linear phase FIR filters one can come up with?

It can be shown that in order for an FIR filter to have linear phase, its impulse response must be symmetric or anti-symmetric. The impulse response length $L$ can be either odd or even. Those two variables lead to 4 combinations, hence the 4 types of linear-phase filters.


Refinement equation and orthogonal wavelets basis

This follows from the orthonormality of shifted versions of $\phi(x)$:

$$\langle \phi(x),\phi(x+n)\rangle =\delta[n]\tag{1}$$

where $\delta[n]$ is the unit impulse. Taking the discrete-time Fourier transform (DTFT) of $(1)$ gives the equation in your question.

Note that the left-hand side of $(1)$ is the auto-correlation of $\phi(x)$ evaluated at integers $n$. The Fourier transform of the auto-correlation is $\left|\hat{\phi}(\omega)\right|^2$, where $\hat{\phi}(\omega)$ is the Fourier transform of $\phi(x)$. Sampling the auto-correlation results in the infinite sum of shifted spectra.

EDIT: The equation can be derived as follows: define the auto-correlation of $\phi(t)$ as

$$p(\tau)=\int_{-\infty}^{\infty}\phi(t)\phi(t+\tau)dt\tag{2}$$

Realizing that $p(\tau)=(\phi_{-} \star \phi)(\tau)$ (where $\phi_{-}(t)=\phi(-t)$), the Fourier transform of $p(\tau)$ is

$$P(\omega)=\left|\hat{\phi}(\omega)\right|^2\tag{3}$$

Sampling $p(\tau)$ at integer $n$ results in a periodic continuation of the spectrum:

$$\textrm{DTFT}\{p(n)\}=\sum_{k=-\infty}^{\infty}P(\omega+2\pi k)\tag{4}$$

With $p(n)=\langle \phi(x),\phi(x+n)\rangle$ and with $(3)$ we finally get

$$\textrm{DTFT}\left\{\langle \phi(x),\phi(x+n)\rangle\right\}=\sum_{k=-\infty}^{\infty}\left|\hat{\phi}(\omega+2\pi k)\right|^2\tag{5}$$


Understanding the relationship between the convolutional feature map and the region proposals in Fast R-CNN

Notation

Lets say you have an input image I with size WxHx3
ROI proposal x0,y0,w0,h0.
When you forward the image through the feature extractor you get a feature map F with size WfxHfXC size.

The relation is a spatial relation
The ROI doesn't affect the number of output channels from the feature map.
The relation is a spatial relation which maps an input ROI to an equivalent spatial patch on the feature maps which has the same number of channels has the feature map. So the ROI in the output feature map will be of size UxVxC the same number of channels the feature map has.
Tracking an input pixel throgh a CNN
Now lets track a the pixel x0,y0 between successive layers of a CNN.
For simplicity i will consider only the following layers types (The same has in VGG 16 which is the feature extractor used in the article, for other layers the same logic can be applied with small variations):

Convolution layer with "SAME" padding, kernel = 3, stride = 1 and denoted by Conv
Pooling layer with some stride = 2 and denoted by Pool

Now what happens when we insert an image to Convolution layer? In other words to which pixel an input pixel x0,y0 is mapped to in the output of the Convolution?
The output size of a convolutional layer is calculated using the following formula: Output size calculation so in our case we will get that the output size is equal to the input size and that the pixel in the convolution response for pixel x0,y0 in the input is mapped to pixel x0,y0 in the output of the convolution.
For a polling layer with stride 2 we will get that if the input size is WxH the output size will be W/2xH/2 so pixel x0,y0 in the input is mapped to floor(x0/2), floor(y0/2)
Summery of tracking (Given our simplified case)

Convolution x0,y0 in the input will be mapped to x0,y0 in the output
Pooling x0,y0 will be mapped to x0/2,y0/2

Converting x0,y0 in image to x0^,y0^ in the output feature map
So if we have a CNN with Conv->Conv->Pool->Conv->Conv->Pool ....
The pixel x0,y0 in the input will be mapped to the pixel floor(x0/2^(num of pool layers)),y0/2^(num of pool layers)
An ROI mapping
Lets take the case in which the feature extractor is VGG16, which the network used in the article.
All the Convolution layers give output size is the same has the input size and the polling size is half the input size.
The input size is 224X224 and the feature map is 7X7 -> the pixel x0,y0 in the input image is mapped to x0/32,y0/32 in the feature map.
The roi x0,y0,w,h is mapped to x0/32,y0/32,h/32,w/32


Audio template matching - power correlation?

If you're cross correlating short time power spectra and disregard phase, then you're not really disregarding phase.

Phase of the global Fourier transform encodes the temporal structure, including the position of the transients, evolution of tones or the local incoherence of noise.

By capturing the time dependence in the moving window of the time-frequency power spectrum, the really important aspects of phase are represented in the temporal evolution of the frequency power density.

So what you are discarding is merely "local phase", which contains a lot of information that is however modified by even the most subtle processes like sound propagation, speaker reproduction, microphone recording, etc. These modifications don't affect the qualitative content of the sound a lot, and a robust sound recognition algorithm should be mostly insensitive to them.

That means discarding the local phase will make your recognition or correlation algorithm more robust and avoid misclassifications due to small inaudible phase errors, while at the same time preserving the total temporal structure of the signal.


I’m guessing your book is on Kalman filters and $u(k)$ is a term in a state variable equation of a linear system.  If so, it is a vector and the components of that vector at some time index $k$ can be correlated with each other but are independent with respect to $k$.

it isn’t a dot product, it is an outer product.


A confusion may arise, for causal systems, from mistaking "having negative signal values (amplitudes)" and "depending on  negative time indices".

A strict memory-less system does depend neither on past (for the causal case) nor on future values, only the current ones, to determine the current value of the output. The output at $n$ only depends on inputs at $n$. For a strict memory-less and causal system, you could not depend on former (negative, in a way relative to the current index) time indices.

By extension, some allow memory-less systems to allow a constant delay, i.e. to hold only one value at a given lag: the output at $n$ only depends on inputs at $n-k$  of a fixed integer $k$: So, globally, memory-less system can afford the following shape, with $k$ any integer:

$$y[n] = f(x[n-k])\,.$$

For realizable causal systems, that cannot look ahead in the future, on:y $k\ge 0 $ are allowed.


If you have reason to believe your signal is sparse in the frequency domain, they you should be denoising in that: use FFTs, not wavelets. Look at the magnitude spectrum and see if it's very spiky, and that will be your cue. If so, attenuate the frequencies with a small response.


How should we read FIR filter phase response graph?

First of all you see that the phase is a piecewise linear function, so it's a linear phase FIR filter. There's a phase jump at half the Nyquist frequency, which shows that the filter has a zero at that frequency. Note that the phase jumps by $\pi$ corresponding to sign inversion.

You can also see what type of linear phase FIR filter it is. Since the phase is zero at $\omega=0$ it must be either type I or type II (i.e., even symmetry). If it were type III or type IV (odd symmetry), the phase would jump at $\omega=0$, and the phase at $\omega=0^+$ would equal either $\pi/2$ or $-\pi/2$. You can distinguish type I and type II by looking at the phase at the Nyquist frequency. A type I filter (odd filter length) must have a phase of either zero or $\pi$ at Nyquist, whereas the phase of a type II filter jumps at Nyquist because a type II filter always has a zero at Nyquist. You can't see the phase jump here (because the figure doesn't go beyond Nyquist) but you can see that the phase equals $-\pi/2$ just before Nyquist. A type II filter's phase always equals $\pm\pi/2$ at $\omega=\pi^-$, and it jumps by $\pi$ at $\omega=\pi$.


Peaking filter (loudness) cutoff frequency in an equalizer

The purpose of the loudness filter is to compensate for the level dependency of the a human's "frequency response". Most music is mixed at pretty high levels, say at 80-90 dB SPL. If you play this back at much lower levels, say 40 dB SPL, it will sound very bass deficient. The loudness filter compensates for this by boosting the bass at lower levels.
In a "good" implementation the loudness adjusts automatically based playback volume setting and content energy: the lower the playback volume the more bass boost is needed to maintain the same perceived spectral balance between low and mid frequencies.
The shape of the loudness filter can be derived from the Equal Loudness contour of humans (see https://en.wikipedia.org/wiki/Equal-loudness_contour). In essence the loudness filter is the difference between the equal loudness curve at the playback level and the equal loudness curve at some reference level, all normalized to 0 dB at 1 kHz (or something like this).
The target filters have a continuously rising bass response. That's not easy to implement with "standard" EQ blocks. A very low frequency peak with the right choice of center frequency and Q can be a fairly good match in the region you care about. You use a 20 Hz peak filter because it has the right shape from 40Hz-100Hz and you don't care much what happens at 20 Hz since there almost no content there.


Several questions about power of periodic signal

You are missing the effect of the integration interval on the integrand function...

If you choose a period; T = 2 , and a periodic signal whose base period is $s_0(t) = t$  in the interval [0,2]... Then the integrand in the shifted interval [-1,1] will be different as given by

$$
s_1(t) = \begin{cases}{ t + 2  ~~~, -1 < t < 0 \\ t   ~~~~~~~~~~ , ~   0< t < 1}\end{cases}
$$

Then you should write the integrals as follows: :
$$
P = \frac{1}{2}\int_{-1}^{1}|s_1(t)|^2dt = \frac{1}{2}\int_{0}^{2}|s_0(t)|^2dt
$$

You should get the same result now...


Performance of MATLAB's `conv2()` vs. `imfilter()`

First, I think conv2 and imfilter will give you the same result if you change the filter option of imfilter to conv instead of corr. imfilter uses correlation to filter images by default that starts from one side of the image, whereas covolution starts from the other, so there may be some small differences in the filter output.

Secondly, my test shows imfilter on 2D image is faster if the filter size if the same as images.

clc
data = (rand(200));
kernel = ones(200)/4e4;
tic
blurredData = conv2(data, kernel,'same');
conv2t=toc
tic
blurredData = imfilter(data, kernel,'same','conv');
imfiltime=toc

>> conv2t =

       2.7425

  imfiltime =

       0.1552

Both imfilter and conv2 are hardware optimized and multi-threaded. For floating point inputs, the performance of them are supposed to be roughly the same. But in most cases in image processing, the filter kernel is linearly separable, imfilter does separation of the kernel into row and column vectors. As a result, they are much faster.  Yet if the filter kernel is (1) not linearly separable, or (2) a small kernel such as a 5*5 filter kernel. conv2 wins in the processing time. I guess the reason is some time is spent on checking the separability of the kernel in imfilter. I am not sure whether there are other different mechanisms that leads to the processing time difference between the two.

Finally, if the filter size is the same as image. fft is always the best one (O(nlogn)):

clc
data = (rand(200));
kernel = ones(200)/4e4;
tic
blurredData = conv2(data, kernel,'same');
conv2t=toc
tic
blurredData = imfilter(data, kernel,'same','conv');
imfiltime=toc
tic
D=fft2(data);
K=fft2(kernel);
blurdata=ifft2(D.*K);
fftime=toc

>>conv2t =

          2.7429

  imfiltime =

          0.1327

  fftime =

          0.0164

EDIT

@geometrikal's detailed analysis:
The reason fft is faster than imfilter with same size images is that imfilter pads the inputs (even with 'same'). For small filter kernels, not much padding is needed for imfilter, whereas to use the fft the kernel will have to be padded to the size of the image. In these cases imfilter is faster.


Sampling is the process of making the x-axis (time) discrete and quantization is the process of making the y-axis (magnitude) discrete.  You can sample without quantization (such as done with an analog sample and hold circuit). Quantization is introduced through rounding or truncation when the sampled analog signal is mapped to a digital representation.

Sampling only causes aliasing distortion to the signal; if there is no energy above the half the sampling rate for real signals with real sampling, or above the sampling rate for complex signals with complex sampling, absolutely no distortion occurs assuming a perfect sampling clock with no hitter/phase noise (both aspects are not physically possible to totally eliminate but can each be reduced to below our level of concern).

Many A/D's, such as successive approximation converters, are implemented with a fast analog sample and hold so that the sampled signal is held while the conversion to digital takes place.

Quantization noise adds noise to the resulting signal, well approximated in most conditions as a uniform white noise: the noise distribution in magnitude is uniformly distributed over the magnitude of a quantization level and is white in frequency to the extent that each sample of the noise is independent of the next.

From that approximation of the noise we can arrive at the SNR estimate for a full scale sine wave to be 6.02 dB/bit + 1.76 dB, which we can increase further by “over sampling”; since the noise added is white, we can filter out the noise in the band we don’t need.

I have detailed that process further at this post: What are advantages of having higher sampling rate of a signal?


The effect of upsampling on DFT coefficients

Note that if you upsample a sequence $x[n]$ by a factor $M$ you have

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

The DFT of $x[n]$ is

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j\frac{2\pi}{N}nk}\tag{2}$$

The DFT length of $y[n]$ is $MN$, so the DFT of $y[n]$ is

$$Y[k]=\sum_{n=0}^{MN-1}y[n]e^{-j\frac{2\pi}{MN}nk}\tag{3}$$

Since only every $M^{th}$ sample of $y[n]$ is non-zero, $(3)$ can be written as

$$Y[k]=\sum_{n=0}^{N-1}y[Mn]e^{-j\frac{2\pi}{MN}Mnk}=\sum_{n=0}^{N-1}y[Mn]e^{-j\frac{2\pi}{N}nk}\tag{4}$$

With $(1)$ we get

$$Y[k]=X[k]\tag{5}$$

Note that from $(2)$ $X[k]$ is $N$-periodic, so the length $MN$ DFT of $y[n]$ is just an $M$-fold repetition of the DFT coefficients of $x[n]$.

A simple example:

x = [1,1];
X = fft(x)
ans =

   2   0

y = [1,0,1,0];
fft(y)
ans =

   2   0   2   0

y=[1,0,0,1,0,0];
fft(y)
ans =

   2   0   2   0   2   0


Can the FFT representation of a signal be modified w/o creating glitches?

The way you do it is you use a window that tapers to 0 at the ends, both before and after the FFT (analysis and re-synthesis windows), along with 50% overlap (needs to be a COLA window), to fade out the glitches and keep the signal continuous.


Advice on filtering a CW signal for processing

The recommendation to use 4* WPM just means that you should estimate the upper limit of the WPM that you expect to receive.  High speed CW starts around 40 WPM (usually too fast to write or type the decoded signal.) So, 50WPM would be a good upper limit which means you could use a bandwidth of 200Hz.

As an alternative to what you are doing, you could use the goertzel algorithm.  That algorithm works like a narrow bandpass with a signal level detector.  You run your signal through it, and it delivers a number that tells you how strong the signal is at the center frequency that you gave the goertzel method.

This web site  gives a clear example of how to write a goertzel method to do what you are trying to accomplish.  Use the tandem goertzel to work in real time.  The example is in C, but it should be clear enough that you can translate it to python.

Varying the RESETSAMPLES value will allow you to change how narrow the filter is.  The higher RESETSAMPLES is, the narrower the filter will be.  The article suggests 200 samples, so I would start there.  If your signal is too noisy and you get a lot of bad morse code symbols, then you can make the number larger.  What I usually do is make it some multiple of the wave length of the frequency I am trying to detect in samples.  So, to detect 1kHz using 44100Hz sampling, I would first try 44 (44100/1000.) If that is too "loose" then 88, then 132 and so forth until I get a filter that is narrow enough.  If you go too high, though, you will start "smearing" the morse code symbols and won't be able to tell when one ends and the next begins.


What are the fundamental differences between speech and music in frequency domain?

One way is to look at the energy below 80 Hz. Unless you have a recording of James Earl Jones, there will be almost none for speech and any music with bass and/or drums is going to have a lot down there.
Note that there is huge variability in music recording and a lot in speech. Creating an algorithm that can reliably distinguish between the two for the majority all recordings is a fairly complicated. Looking at the "average spectrum" isn't that useful since the any individual instance can (and does) look completely different than the "average spectrum". This is actually a good candidate for machine learning although many deterministic algorithms do exist as well.


Algorithm to compute the Variance of the signal?

You want the sample variance:

$$
s^2 = \frac{1}{N-1} \sum_{n=0}^{N-1} (x_n - \mu_x)^2
$$

where the mean, $\mu_x$, is given by

$$
\mu_x = \frac{1}{N} \sum_{n=0}^{N-1} x_n.
$$

The $\frac{1}{N-1}$ factor should be used to account for the bias that is otherwise there if $\frac{1}{N}$ is used. One way to think about this is the sample variance is computed from $N-1$ differences between samples, rather than $N$ samples.


The cyclic prefix is to make sure that any multi-path interference (or other process similar to a linear time-invariant filter in the transmission channel) acts as a circular convolution on the FFT data frame, thus not affecting the orthogonality of the data channels within the FFT bins.  It also allow some slop in the receiver's symbol clock.

Thus the cyclic prefix has to be a circular repetition of the IFFT data, either before or after the IFFT block, and long enough to be as long or longer than the maximum expected multi-path delay time differential in the channel, plus the maximum of any timing error between the received data block and the receiver symbol clock.  The N/4 is just a default estimate of this time for some common channels, but longer may be required.

Anything other that a circular repetition can cause multi-path delays to act as a non-linear spectral spreading function (relative to the FFT frame) which will affect channel orthogonality, and thus S/N related reliability.

Repetition (circular or otherwise) also helps to statistically even out the PAP (peak to average power) ratio over the extended frame.  And prefix repetition may help push any transmitter anomalies due to sudden changes in this ratio forward of the FFT frame.


has some compensation in the signal chain to compensate for the DC spike which would appear at 370MHz the center frequency.

Yes, a typical problem of IQ demodulators is the LO leakage/DC spike. You really can't do much about that. You don't mention the SDR you're using, but many can be used with offset tuning, completely moving the LO out of your band of interest and thus completely solving that problem

  But unfortunately this has caused a dip at the centre frequency.

Why is that a problem? Looking at your spectra, I'd clearly say you've got four subcarriers working there, and the centre frequency happens to end in the gap between the center two. There's really not much you'd need to fix.

  Now since the Tx channel is varying the signal power fluctuates which is shown in image 2 and image 3. The null at 370MHz is at about -22dB but this null should be at about -17dB with the other nulls.

Why? It's a null. There's no data here. Simply leave it alone!

  I want to be able to dynamically compensate for the 370MHz null during varying channel conditions bringing it up to -17db

That's a bad idea. You'll break more than you fix when trying to equalize based on a cancelled component. You get overproportional noise amplification that way!

Instead, you'd want to base your equalization on something that comes through well.

This is exactly the discussion you probably read when you first were introduced to frequency-domain equalization: Zero-Forcing looks easy, but will amplify the noise on channel-suppressed frequencies. You'd be doing even worse: you'd be selectively taking the noise on a suppressed frequency and use that to equalize. That means you'll multiply a lot of noise onto your whole spectrum. Bad, bad idea.

  also compensate for the slope on the flat tops adjacent to the 370MHz null. I have been trying of fit a line to the top of the flat tops with little success.

Why should there even be a flat top here, or a straight line? This is clearly a multipath channel, and a curve fit isn't what you'd use to equalize that.

From the look of your signal, I'd say this is probably one larger or four low-subcarrier count OFDM transmissions; if that is the case:

The whole point of OFDM equalization is done by the OFDM receiver on each subcarrier individually, and is not something that you should even do before.

All in all, the frequency selectivity of this channel looks pretty benign.


First of all, the title of your question is misleading because onset detection is very different from envelope detection. The onset detection function won't tell you how long a note is played or how slow its attack is... the envelope will. If the problem you are trying to solve is envelope detection, you might get inspiration from onset detection literature, but at some points the methods diverge. Onset detection extracts a "cartoon" version of the envelope in which all sudden changes are emphasized to make their detection easy.

  I think we then alter those coefficients in such a way that the new
  Fourier transform coefficients are now balanced perceptually(ie the
  coefficients are now what a person's ear would perceive them to be)

This is done by summing the energy (squared magnitude) of the FFT across 40 ranges of frequencies defined by the mel scale - usually with a triangular shape to emphasize a different center frequency.

Check point 4 of my answer to this question.

  Why and how would we calculate the first order difference?

This is the $1 - z^-1$ filter.

So the first order difference is simply the energy in a given Mel band at time t minus the energy in the same Mel band at time t-1.

We are interested in detecting onsets. Onsets are mostly characterized by a sudden rise of energy across the entire spectrum. A positive value of the first order difference indicate that energy is rising in one band. When summing across all bands, this gives a value that is very likely to be high when there is a note onset.

We do this detection across 40 frequency bands rather than in the time domain to detect events that would be masked by loud sounds in the time domain. For example, imagine a recording in which a high-hat cymbal sound is played on top of a loud, droning bass synth. In the time domain, the envelope of the signal will be a solid line with just a tiny wiggle when the high-hat plays, because the most dominant sound is the bass. But in the frequency domain, the bass note will only occupy the lowest frequencies, and it'll be easy to notice the attack of the high-hat sound, extending from mids to high frequencies.

  Why do we need the signal to be locally zero-mean?

It makes it easier to detect onsets by thresholding. The high-pass filter (or subtraction of median, which is another commonly used detrending method) will make sharp peaks in the detection function more prominent.


Laplace transform of product of signal and impulse train

$\delta(t)$ is indeed the Dirac delta impulse, which is not defined by its values (because it is no ordinary function), but by its properties under an integral.

The first expression in your question is a standard model for sampling a continuous-time signal, and since it's a mathematical model it does not represent actual physical sampling. You can replace the Dirac delta impulse by other impulse-like functions. If you use a rectangular impulse you get a zero-order hold. However, if you're not interested in any effects of non-ideal sampling, multiplication with a Dirac comb is a convenient and useful model for sampling.


How to get more accurate results with OpenCV

This termination criteria tells the algorithm to stop when it has either done
20 iterations or when epsilon is greater than .3.  You can play with these parameters for speed vs accuracy but these values work pretty well in many situations.

So if you want to improve accuracy then you will get more accurate value with more iterations but the time taken for doing such iterations will be improved , I mean speed will be reduced.


How to estimate the signal-to-noise ratio of a waveform?

If the spectrum of the signal and the spectrum of the noise do not overlap, you might be able to integrate or sum the energy in each of the two frequency ranges, and take the ratio.


Every codeword for the channel encoding must have the same length.

That is not inherently true, but it's usually handy to consider the information that goes into the channel coder to be of constant length, especially when considering block coders.

  Why is this? If the source from which the channel is reciving its input used a Huffman code, for example, lengths might be different for each symbol. Isn't that a contradiction?

No. The point is that channel and source coding are two independent things. You take source data, and the source coder maximizes the average information per bit (by converting input words of constant length into source-coded words of varying length). Then the source-coded words are again simply considered as sequences of bits (or to be more general, sequences of fixed-length units of information) and channel-coded, optimizing robustness.


PyWavelets CWT: resampling vs recomputing wavelet

The normalization is indeed by 1 / sqrt(scale), and it's an L2-norm; the trick's in the scale wavelet.
I'll use wavelet='morl' throughout. Prior to integrating, we can inspect the wavelet here; it's returned by wavelet.wavefun, which is binary-compiled, but after some guesswork, I found to it to match exactly with
scipy.signal.morlet2(1024, 1024 / 16) * sqrt(1024 / 16) * np.pi**(.25)

which is, from source, using Wiki's notation, $\psi(t) = \psi_{\sigma}(t/a)$, where $a$ is scale, and
$$
\psi_{\sigma}(t) = e^{j\sigma t} e^{-t^2/2} \tag{1} \label{1}
$$
(scale and $\pi^{-1/4}$ cancel). This is what's integrated via cumsum(psi) * step, and then resampled for all scales.

Resampled vs. recomputed
What exactly is this resampling doing in in terms of Eq 1? Is it just a higher resolution of the wavelet at the same scale, or is it equivalent to recomputing Eq 1 at each scale? Conveniently enough, latter, but only approximately, and approximation degrades substantially for small scale (-- code1):

Notice from code1, however, the recomputed wavelet:
Ns = len(int_psi_scale)
w = real(morlet2(Ns, Ns / 16) * sqrt(Ns / 16) * np.pi**(.25))  # repeats first blob
w /= scale

The recomputing includes 1 / scale, which, together with * sqrt(scale), equals 1 / sqrt(scale). Mystery solved.

Your code is wrong, where's * step?
Replaced by 1 / scale. How?
In MAE code, notice that for scale=64, we have int_psi_scale == int_psi, which == cumsum(psi) * step. For w_int, we do cumsum(w) / scale. And 1 / scale is... == step. Thus, the pre-integrated wavelet, psi, is just w at scale=64 (in morlet2 code above, 1024 / 16 == 64, checks out), and step happens to be ... conveniently? == 1 / scale when integrating.
Then, why is 1 / scale there? Unclear. Two possibilities in mind: (1) preserving the wavelet's norm at integration; (2) scaling the wavelet, independent of integration.

If the wavelet was L1 or L2 normalized before integration, either will be preserved. This is from chain rule; just replace $f$ with $\psi$, and $k$ with $1/a$:

$$
\int f(k x) dx = \frac{1}{|k|} \int f(x) dx
$$

This seems likelier, as the diff later is closely equivalent to undoing the integration, defeating the purpose of (1). Why rescale the wavelet? Normalization - see here.

You cheated earlier; there is no w /= scale
True, the code actually shows w_int = cumsum(w) / scale, but the two are exactly the same. It's thus the earlier question of where 1 / scale "belongs", or "comes from". This is answered here, and in other part below.

Why step == 1 / scale at integration? (-- for reference, from here (in code1, $n$ is x)):

Just a coincidence, or is step, along $n_i$, carefully crafted to yield the convenient resampling properties, which in turn demand step = 1 / scale? Don't know, might update answer later.


Why using the impulse respond of an LTI to characterize its output?

For the two special classes of systems:
$$ \textbf{Linear time and shift invariant systems } $$ An LTI and LSI system is one that obeys superposition and time and shift invariance. The LTI systems are described by linear ordinary differential equations and are your analog systems. The LSI systems are described by your linear constant coefficient difference equations. There are many analogues between them, as is to be expected as the discrete systems can be taken to be the discrete versions of the analog ones.

The dirac pulse and unit impulse sequence completely characterise the systems in both the time and frequency domain.

Let us first define what a system is, it is a relation that takes in an input signal $ x(t) $ or $ x[n] $ and produces an output $y(t)$ or $y[n]$. As is the case for analogue and discrete ones. The system is thus a relation, or an operator on the input to produce the output.

$$
y(t) = \textbf{T} \big \{x(t) \big \}
$$

$$
y[n] = \textbf{T} \big \{x[n] \big \}
$$
We define the response to a dirac pulse and a unit impulse respectively for system classes 1 and 2 as:
$$
h(t) = \textbf{T} \big \{ \delta(t) \big \}
$$

And the unit impulse response for discrete linear shift invariant systems as:
$$
h[n] = \textbf{T} \big \{ \delta[n] \big \}
$$

The impulse response for these two classes of systems is hence, the output of such systems when the input is a dirac pulse or a unit impulse sequence. Note the following properties of the dirac pulse and unit impulse sequence:

Every function can be expressed as a superposition or sum of scaled and shifted dirac pulses or unit impulses.
$$
y(t) = \displaystyle \int_{- \infty}^{\infty} y(\tau) \cdot \delta(\tau - t) \,\,\,\, \text{d}t
$$

$$
y[n]= \displaystyle \sum_{ k \to - \infty}^{k \to \infty} x[k] \cdot \delta[k - n]
$$

The continous time fourier transform and the discrete fourier transform of the dirac pulse and unit impulse are both 1, hence, in the frequency domain, the pulse subjects your system to all possible frequencies, giving a fingerprint of the system.

$$
\textbf{The transfer function of a system (frequency domain or s domain) }
$$

It is the fourier transform of the impulse response, or the laplace transform of the impulse response (s domain).

$$
H(j \omega) = \displaystyle \int_{- \infty}^{\infty} h(t) \cdot e^{-j(\omega t) } \,\,\,\, \text{d}t
$$

And

$$
 H(s) = \displaystyle \int_{- \infty}^{\infty} h(t) \cdot e^{-s (t) } \,\,\,\, \text{d}t
$$

Now we have to talk about convolution and multiplication time duality.
Multiplying two functions in the frequency domain is equal to convolving them in the frequency domain.
Convolving in the time domain is the same as multiplying in the frequency domain.

Now to put it all togther, since your input function can be expressed as a sum of scaled and shifted dirac pulse or unit impulses, your output, due to linearity and time invariance, would be equally scaled and shifted impulse responses:
$$ a h(t - t_{0} ) = \textbf{T} \big \{ a \delta(t - t_{0} ) \big \} $$
This is where convolution comes in, you express the input function as a sum of scaled and shifted dirac pulses or unit impulses, you then input these into your system, and to get the output you sum up or integrate all the scaled and shifted impulse responses, to get your output.

$$
y(t) = \displaystyle \int_{-\infty}^{\infty} f( \tau) \cdot h(t - \tau) \,\,\,\, \text{d} \tau
$$

$$
y[n] = \displaystyle \sum_{k \to - \infty}^{ k \to \infty} x[k] \cdot h[n - k]
$$
And in the frequency domain due to the theorems this is equal to
$$
Y(s) = H(s) \cdot X(s)
$$


Quite often, performing a time-delay estimation in a transformed domain can help,as long as the signal gets better concentrated, and/or the noise is more spread.

So performing the CC or GCC in a time-frequency or time-scale (wavelet) domain could be useful, as long as you can select the frequency  band or the scale in which the correlation will be meaningful.

If the noise characterics are unknown, you can benefit from wavelet whitening capabilities.

You can check for instance:

On Wavelet Denoising and its Applications to Time Delay Estimation, 1999

For Matlab codes:

Wavelet Cross-Correlation for Lead-Lag Analysis
modwtxcorr: Wavelet cross-correlation sequence estimates using the maximal overlap discrete wavelet transform (MODWT)


What is the algorithm to do frequency modulation?

What your plot shows is noise from varying size discontinuities in the waveform.

When you change a modulating frequency, you have to make sure the resulting FM waveform does not have discontinuities.  You have to do this by incrementally changing the phase of the argument to the sin() function, not by naively multiplying the modulating frequency by time, which will jump the phase by some (psuedo)random angle when the frequency changes.


Are all LTI systems invertible? If not, what is a good counterexample?

You need to define what you mean by "invertible". Do you mean invertible by a causal and stable system? If yes, then any system that is not minimum-phase is not invertible (because the inverse system can't be causal and stable).

Example of a system that cannot be inverted by a causal and stable system: a simple delay $y(t)=x(t-T)$, $T>0$, could only be inverted by a non-causal system.

If by "invertible" you mean a system that can be inverted by a system that is possibly not causal and/or not stable, then it's still straightforward to find systems that can't be inverted according to this broader criterion: just take a system with one or more zeros in its frequency response. The information about the input signal at the frequencies where the system's frequency response is zero cannot be recovered by any other system.

Example of a system that cannot be inverted by any system: an ideal low-pass filter cannot be inverted by any system, because any information above the filter's cut-off frequency is lost and cannot be recovered.

EDIT: the latter definition of invertibility (where causality and stability of the inverse system are ignored) is just the definition of an injective mapping. In terms of systems this means that no two distinct input signals generate the same output signal. This category is what is also described in Laurent Duval's answer. Note, however, that in practice we usually prefer inverse systems that can also be realized, i.e., we require stability and causality, and then the first definition of invertibility given above is appropriate.


Confusion about hough transformation for ellipse

It's a little mis-leading and confusing.

The usual way to describe an ellipse is using cartesian ($x,y$) coordinates.  Another way to represent an ellipse is using polar coordinates ($r,\theta$) (where $r$ is the distance (radius) from the origin, and $\theta$ is the angle.

What the equation is saying is that the cartesian coordinates of the ellipse are parameterized by $\theta$ in the two equations. That is, as you change $\theta$ from 0 to $2\pi$, the values for $x$ and $y$ will lie on the ellipse, for constants ($a_0,b_0,a_x,b_x,a_y,b_y$).


For the sake of putting some numbers to this question, I implemented a basic histogram of gradient from scikit-image (skimage.feature.hog). Here is the timing data for HOG with default parameters applied to the skimage.data.astronaut image in b&w and rescaled to have the given dimensions:

Image dimensions......(102, 102)
3.49 ms ± 26.3 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Image dimensions......(256, 256)
27.5 ms ± 295 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
Image dimensions......(512, 512)
123 ms ± 2.88 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
Image dimensions......(768, 768)
284 ms ± 9.18 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
Image dimensions......(1024, 1024)
509 ms ± 8.62 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
Image dimensions......(2048, 2048)
2 s ± 15.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

And here is timing data for the inverse Fourier transform (np.fft.ifft2) on the Fourier transform of the same set of images. (Note that taking the Fourier transform is not included in the timing results)

Image dimensions......(102, 102)
684 µs ± 32.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
Image dimensions......(256, 256)
2.51 ms ± 26.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Image dimensions......(512, 512)
21 ms ± 519 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Image dimensions......(768, 768)
42.5 ms ± 621 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
Image dimensions......(1024, 1024)
85.9 ms ± 1.47 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
Image dimensions......(2048, 2048)
381 ms ± 9.18 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Graph of the timing data


How are Object Recognition and Content-based Image Retrieval interlinked?

The difference is not so much in the algorithms used, but in the objectives.  Object recognition is a problem of naming an object or objects depicted in an image.  Content-based Image Retrieval is the problem of finding images in a database which match a user's query, which may be in the form of an image or in form of text.

I would say that CBIR is a more general problem than object recognition.  If anything, it is CBIR that builds on object recognition, and not the other way around.  If you have an image retrieval system, and the query is an image of a cat, then your system would need to recognize that this is an image of a cat (object recognition), and then find other images of cats in its database.  Presumably to be able to do that it would have to have performed object recognition on all of its images in advance, and tagged the ones that depict cats accordingly.

Note, that you can also have a query that may not require object recognition.  For example you may ask the system for images of urban landscapes, or images of indoor scenes.  This is a problem of "scene recognition" and it is typically approached  by using  some sort of a global image descriptor (e. g. scene gist), rather than trying to find all the object in the image.

Also note that an object recognition system must give you a single answer, i. e. the name of the object in the given image, and that answer is either right or wrong.  On the other hand, an image retrieval system gives you multiple possible matches for your query, and you would consider it successful even if only some of the top matches are good.  So there is a big difference in how you might evaluate an object recognition system vs. an image retrieval system.


Absolute summability implies that the Fourier transform converges uniformly, so this means that the unit circle is contained in the region of convergence (ROC).

Now, the pole at 0.5 is the limit of the ROC, so the ROC is from 0.5 outwards (in order to include the unit circle). This means that the signal is right-sided.

Finite support signals cannot have poles (other than $z=0$ or $z=\infty$).


From the code in Olli's comment, the phase is going to be

beta = t1 / log(f1 / f0)
phase = 2 * pi * beta * f0 * (pow(f1 / f0, t / t1) - 1.0)

So to ensure you have a zero crossing at the end of your signal, you need to choose a combination of t's end value and t1 for which the phase is modulo $2\pi$.

From your code, the last value of t will be 44099/44100. So the beta value is:

beta = 1 / log(1000)

and the phase of the last sample in the signal will be

phase =  2 * pi * 1/log(1000) * 20 * (pow(1000, 44099/44100) - 1.0)

which will be close to modulo $2\pi$ if you have a ratio of f1 to f0 equal to 1000.007748 instead of 1000 (after a bit of playing in Excel).  Try setting f1 to 20000.15495 instead of 20000.

EDIT : I made a mistake in the above and have corrected it. Here's a Google Doc that has the spreadsheet I've used (read only) for the correction.

Screenshot highlighting the close values below.


System Identification with a Limited Bandwidth Input Signal and Region of Interest

The reason x[n] must be white is because the solution will effectively spectrally weight the channel response based on the amount of energy present in each spectral frequency location. A white noise source provides equal weight to all frequencies. If energy is not present in any particular frequency bin, a proper solution cannot be found for that frequency.

If you only care about a small band of signals, then I would argue that you can still use the least squared approach. The reason for this is to consider a system that provides band pass filtering: for such a system I could excite the input with a white noise source and use the least squares approach (Wiener-Hopf equations) of comparing the input and output signals which would accurately provide a least squares estimate of the channel. So if you have a band limited signal, as long as that signal is white over your band of interest, it will still provide the accurate solution over that band (and you ignore everything else).


How to measure the agreement between to curves?

If the signals are not aligned, yet you have a clue that they more or less "mean" the same thing, or refer to similar data, you could use dynamic time warping (DTW) algorithm to obtain a better correspondence (then simply taking the value at the same location). You could measure RMS, MSE or whatever you like, using those correspondences. For DTW, you might want to check:
http://en.wikipedia.org/wiki/Dynamic_time_warping

A good way to elaborate this idea is utilized in Earth Mover Distance (EMD), which computes the minimal effort to bring the signals to alignment, as a distance measure. EMD is presented here:
http://ai.stanford.edu/~rubner/emd/default.htm

EMD gives you a direct distance, which you might use for further analysis.


The output of your filter is what is sometimes called band-limited white noise. In
your particular case, the autocorrelation function of the output noise is a sinc
function whose zeroes are every $100$ microseconds, that is, samples taken at
the rate of $10^4$ samples/second are uncorrelated.  Your samples at $3\times 10^4$
samples per second are closer together and thus are correlated.  In fact,
the correlation coefficient between successive samples is
$\displaystyle \frac{\sin(\pi/3)}{\pi/3} = \frac{3\sqrt{3}}{2\pi} \approx 0.827$

Note added after question was edited: Samples taken $0.03$ milliseconds apart
are at a frequency of $33.33\ldots$ kHz, not $30$ kHz as you say in the
part of the question that you typed in yourself. Regardless, the answer (B) is
incorrect and the reasoning given in support of answer (B)
is bogus. What is being sampled is
not white noise but filtered white noise, and Answer A is correct.

For (B) to be the correct answer,
the sampling rate must be a sub-multiple of $10$ kilosamples per second,
that is the samples must be spaced $100$ microseconds apart
or an integer multiple of $100$ microseconds
apart.  Sampling every $0.3$ milliseconds, that is, every $300$ microseconds
meets this criterion; sampling every $0.03$ milliseconds does not.
Now sampling at $10$ kilosamples/second  makes the samples uncorrelated
(answer (C)) and to get from this to the stronger result that
the samples are statistically independent (answer (B)), we need
further assumptions about the noise. The standard assumption is
that the noise is Gaussian.

  In summary, (B) would be the correct answer
  if the sampling were done every $0.3$ milliseconds (since the problem
  statement already includes the assertion that the noise is Gaussian)
  and (C) would be the correct answer if it did not say that the noise
  is Gaussian since we could not the make the specialization from
  uncorrelated to independent.
  But when the sampling is done every $0.03$ milliseconds, then (A)
  is the correct answer.


I don't know about real time application normally what I am doing is to calculate $\mathrm{PRD}$ first and calculate $\mathrm{SNR}$ from $\mathrm{PRD}$

$$\mathrm{SNR}=-20\log_{10}(0.01\mathrm{PRD})$$

$\mathrm{PRD}$ -  Percentage Root Mean Square Deviation


What's the Difference between Convolution Kernel and Point Spread Function in the Context of Image Convolution?

The PSF (Point Spread Function) is the system response to Impulse Signal (Point).

If your system model is LSI (Linear Spatially Invariant) then the output image of the system is applying the PSF as a convolution kernel on the input image.

Yet, the PSF is just the response of a system to a certain input.
It doesn't describe the whole system unless it is an LSI system.

If you model for the Deconvolution problem (Often called the Degradation Model) is LSI, then what you estimate is indeed the PSF.


Removing specific frequency from a signal

Below is a picture that has:

Top graph has the noiseless 20Hz sine wave sampled at 100Hz (blue) and the noisy sine wave (red).
Bottom graph has the actual noise signal added (red) and the noisy signal filtered with an IIR notch filter with a notch centred at 20Hz (blue).

Apart from the group delay between the original noise and the output of the notch filter, the two are very similar.

R code to implement this is below the graph.

The frequency response of the notch filter is as below.

R Code

#26866

f <- 20
fs <- 100
phi <- 2*pi*0.12987892374

T <- 100
t <- 0:(T-1)

signal <- sin(2*pi*f/fs*t+phi)
noise <-  rnorm(T,0,1)
noisy_signal <- signal + noise

alpha <- 0.9

num <- c(1, -2*cos(2*pi*f/fs), 1)
den <- c(1, -2*alpha*cos(2*pi*f/fs), alpha*alpha)

filtered <- signal::filter(num,den,noisy_signal)

par(mfrow=c(2,1))
plot(noisy_signal, type="l", col="red",  lwd=10)
lines(signal, type="l", col="blue", lwd = 3)
title("Noiseless signal (blue) and noisy signal (red)")

plot(filtered, col="blue", lwd = 4)
lines(noise, col="red",  lwd=3)
title("Filtered signal (blue) and original noise signal (red)")


On the meaning of s-plane and it's link to a transfer function

In the Laplace transform, being on the line $s = 0 \pm j\omega$ means you are looking at the pure sinusoids contained in the signal. This is exactly as you mentioned in your understanding of Fourier analysis. So this means that you are looking at the time-domain terms:

$$e^{j\omega{t}}$$

Which again represent the pure sinusoids composing your signal.

However once you are off this line and have some $s$ such that: $$s = \sigma \pm j\omega$$

Not only are you looking at the $e^{j\omega{t}}$ sinusoids terms, you are now looking at time domain terms:
$$e^{\sigma{t}}e^{j\omega{t}}$$

When talking about a stable system, poles being on the left-hand-plane of the $s$-domain means that $\sigma < 0$ so the sinusoids you are considering contained in the signal now have a damping envelope. If $\sigma$ was positive, the poles would be placed in the right-hand-plane and the signal would explode to infinity as time went on. This is why we want poles to stay out of the right-hand-plane if we want a stable system.

The takeaway is that in Fourier analysis you are considering the pure sinusoids composing your signal. With the Laplace transform you are considering sinusoids with that exponential term attached. In the case for a stable or oscillating system, you are considering the damped sinusoids ($s = -\sigma \pm j\omega$), decaying exponentials ($s = -\sigma$), pure sinusoids ($s = \pm j\omega$), and constant amplitude ($s = 0$) components composing your signal.

This is what makes Fourier analysis appropriate for identifying frequency content and the Laplace transform appropriate for analyzing stability and performance parameters.


Remember that the signal is not the samples; the signal passes through the samples and exists between them, too.  So the maximum value of one chunk cannot be less than the minimum value of the next; they can only be equal, as the signal passes from a lower sample in one chunk to a higher sample in the next chunk.

Plotting a solid line from the min of a chunk to the max sample value of a chunk is not a very accurate representation of the signal anyway (because of intersample peaks and not showing the density of the wave).  If that's good enough of a representation for you, I'd just approximate it by drawing to the midpoint in this case:

find max of chunk 1
find min of chunk 2
if max(chunk 1) < min(chunk 2):
    in column 1: draw line from min of chunk 1 to average(max(chunk 1), min(chunk 2))
    in column 2: draw line from average(max(chunk 1), min(chunk 2)) to max(chunk 2)

If you want to get fancier, you could upsample the signal to get a better idea of the intersample peaks, and then do the same thing based on those.

And the best is to use brightness to show the density of the wave at different amplitudes, but that's significantly more complicated: https://flic.kr/p/7S8oHA


How to add noise to time series data in Python?

Some comments:

I'm not sure if this is the right forum for python questions.
You left some relevant information, like the type of variable your data is stored.
I assume that "la place" is laplace distribution.

Please see if this helps, and make the necessary adjustments:

import numpy as np
def fn_addnoise(data):
    i = len(data)
    # create 1D numpy data:
    npdata = np.asarray(data).reshape((i))
    # add uniform noise:
    u = npdata + np.random.uniform(size=npdata.shape)
    # add laplace noise:
    p = npdata + np.random.laplace(loc=0.0, scale=1.0, size=npdata.shape)

    print npdata
    print u
    print p
    return u,p
# 'data' might be a tuple...
data = (1,2,3,4,5,6)
# or a list...
data = [1,2,3,4,5,6,7]
fn_addnoise(data)


Eigen Value and Eigen Vectors in Image Processing

I will try giving you some intuition.
The SVD says each matrix can be decomposed into 3 operations - Rotation, Stretching (Scaling) and the another Rotation.

What matters is which directions are scaled and how.
Directions are vectors (Pointing some direction).

The SVD has many uses in Linear Algebra.
One its most known use is Low Rank Approximation of a matrix.
Low rank meaning to try represent the same data with lower number of degrees of freedom.

It is immediately rings a bell for compression and dimensionality reduction (Also the SVD is practically the basis for Karhunen Loève Compression).

Image Compression

Think of decomposing the image into patches of $ n \times n $ (Where usually $ n = 8 $ or $ n = 16 $).
Now build a matrix of those patches (Each patch as a column).
Then, using the SVD you will be able to find a low rank representation of this image which is basically a compression.
Why isn't it used?
Because this is the most effective compression in $ {L}_{2} $ sense yet the base it uses is adaptive to the image. Namely with each compressed image you'll have also to deliver the base which hurts the practical compression ratio.
But this is the idea.

Eigen Faces

This is a nice implementation of idea using SVD.
You can read about in the great article on Wikipedia - Eigen Face (Including Code).
Basically, again, you can use the low rank representation either to compress the data base of faces and / or use it as a method to recognize faces.

Compressed Sensing

The SVD is heavily used in the Compressed Sensing frame work.
For instance have a look on the Dictionary Learning method called - K-SVD.

The Wikipedia article about the SVD is really great place start reading about it - Singular Value Decomposition (SVD).


How to Extract Nonorthogonal PCA Principal Components

Let's say x,y and z are your basis vectors and you mix them to get what you are calling your mixed sources. If your basis vectors were (x-y), y, and (z+y) then you could mix them in a different way and end up with the same mixed sources.

Basically, there are many ways to come up with a set of basis vectors for your signals. You need to say what kind of basis you want. It is not enough to say that the basis vectors are non-orthogonal.


Terminologies: Adaptive, recursive and iterative

Q1: In short, is adaptive=recursive=iterative?

Filtering is applying a filter $f$ to an input signal $x$ to get an output signal $y$:

$$
y(t) = f(x(t), t)
$$

The filter $f$ is called iterative if its next value can be calculated from other time-dependent values of either $x$ or $y$ or both.

The filter $f$ is called recursive if it depends on $y$ for some time $s \le t$ (where $t$ is the current time).

If the way $f$ uses past values of $y$ to filter $x$ changes, then the filter $f$ is also adaptive.

For example, the filter:
$$
y(t) = x(t)
$$
is iterative, but not recursive nor adaptive.

For example, the filter:
$$
y(t) = x(t) + \alpha y(t-1)
$$
is both iterative and recursive, but not adaptive.

For example, the filter:
$$
y(t) = x(t) + \alpha(y,x) y(t-1)
$$
is both recursive and adaptive (the filter coefficient $\alpha$ depends on some function of past inputs $x$ and past outputs $y$).

  Q2: Where is EM and Kalman Filtering approaches categorized under and why?

Kalman filtering is definitely iterative and recursive as the new updates will depend on past values of the input and output (well state estimates).

The standard implementation of the Kalman filter is not adaptive.  However, there are simple extensions that make it adaptive.


Block diagram for a complex impulse response

A digital system with complex input and complex output has a block diagram in which the various arrows carry complex signals, each register holds a complex number, each multiplier computes the product of two complex numbers, each adder computes the sum of two complex signals, etc. Consider, for illustration, a simple FIR filter with complex input $x$ and complex output $y$ where
$$y[n] = x[n] + a\cdot x[n-1]\tag{1}$$
with $a$ being a complex constant.  The corresponding block diagram is as shown below:

But, in real life, complex-valued signals actually consist of two separate real-valued signals traveling on two separate wires and it is only in our minds that we think these two real-valued signals as constituting a complex-valued signal (In software or MATLAB, each complex sample is stored as a pair of real numbers (two words)).  So, if we express complex input signal $x$ as $x=x_r+jx_i$ and complex output signal $y$ as $y= y_r + jy_i$ where $x_r, x_i, y_r, y_i$ are real-valued signals, what we get are  two circuits each with separate (real-valued) input and separate real-valued output, and internal cross-connections between the two circuits. We have that
\begin{align}
y[n] &= x[n] + a\cdot x[n-1]\\
&\Downarrow\\
y_r[n] + jy_i[n] &= x_r[n] + jx_i[n]
+ \big(a_r +ja_i\big)\big(x_r[n-1]+jx_i[n-1\big])\\
&\Downarrow\\
y_r[n] &= x_r[n] + a_rx_r[n-1]-a_ix_i[n-1]\tag{2}\\
y_i[n] &= x_i[n] + a_rx_i[n-1] + a_ix_r[n-1]\tag{3}
\end{align}
We see that the circuit with input $x_r$ has output $y_r$ that depends on its inputs $x_r[n]$ and $x_r[n-1]$ (as it should) in exactly the same way as $(1)$ but also includes the cross-connection term $-\alpha_ix_i[n-1]$ which depends on $x_i$, the input to the other circuit. Similarly, the circuit with input $x_i$ has output $y_i$ which depends not just on $x_i[n]$ and $x_i[n-1]$ (as it should) in exactly the same way as $(1)$ but also includes the cross-connection term $\alpha_i x_r[n-1]$ which depends on $x_r$, the input to the other circuit. Thus, the block diagram above can be duplicated with some minor changes to labels (e.g. $x[n]$ is replaced by $x_r[n]$ in one circuit and by $x_i[n]$ in the other circuit, and similarly for $x[n-1]$, $a$ and $y[n]$) and the cross-connections added.  I will leave the actual task of sketching the full diagram to those interested in crossing every i and dotting every t.

And that's it.  No messing around with $\sin()$ and $\cos()$ or orthogonality, no phase-shifting by $\pi/2$, or time-shifting


Dealing with Random Errors in Sampling Rates?

You might read up on "Aperture Jitter". You are correct that higher frequencies are affected more than lower ones. Here is a good reference:

http://www.analog.com/static/imported-files/tutorials/MT-007.pdf


Is Periodic Convolution Same as Circular Convolution?

You did not give specific details on how you define periodic or circular. Assuming standard definitions, the answer is yes.

Standard infinite-support data is treated with standard convolution. When the support  is finite, the "outer unknown samples" require assumptions.

Zero-valued outside samples is a common assertion. Another option, often used with Fourier or wavelet transforms, consists in arranging finite-support data on a ring or torus, the first sample following the last one.

It is often termed periodic, cyclic or circular, with a similar meaning.


Spectral density of modulated noise process

This is a standard problem in (among others) the theory of analog or digital communication systems. If $F(t)$ is real-valued, then the standard gimmick is indeed - as already pointed out in a comment by Dilip Sarwate - to add a random phase to the carrier phase. This is not only a trick but it actually reflects reality in the sense that we often do not know (and do not care to know) the exact phase of the carrier with respect to the signal $F(t)$. So what we actually want is NOT the power spectrum of $F(t)\cos(\omega_0 t)$ (which does not exist), but the power spectrum of the process $F(t)\cos(\omega_0 t+\theta)$ with some unknown phase $\theta$, the obvious model of which is a random variable with a uniform distribution in $[0,\pi)$. This random phase is of course constant over time, but its value is chosen at random reflecting our uncertainty about the origin of the time axis. With this random phase the resulting process is wide-sense stationary (WSS) (if $F(t)$ is), and the power spectrum looks exactly as one would expect.

Now comes the interesting part. If $F(t)$ is a WSS complex-valued process and if we're interested in the power spectrum of the process

$$G(t)=\text{Re}\{F(t)e^{j\omega_0 t}\}\tag{1}$$

then we could also use the random phase trick, but under some reasonable assumptions we don't even need it. The process $G(t)$ in (1) occurs e.g. in systems using quadrature modulation. It is straightforward to show that the autocorrelation of $G(t)$ is

$$\begin{align}R_G(t,\tau)&=E\{G(t+\tau)G^*(t)\}=\\&=\frac12\text{Re}\{R_F(\tau)e^{j\omega_0\tau}\}+\frac12\text{Re}\{R_{FF^*}(\tau)e^{j\omega_0(\tau+2t)}\}\tag{2}\end{align}$$

which becomes independent of $t$ if $R_{FF^*}(\tau)=0$, i.e. if the cross-correlation of $F(t)$ and $F^*(t)$, vanishes. If $F(t)=A(t)+jB(t)$, it can be shown that this is the case if

$$R_A(\tau)=R_B(\tau)\tag{3}$$

and

$$R_{AB}(\tau)=-R_{AB}(-\tau)\tag{4}$$

So if conditions (3) and (4) are satisfied, we don't even need a random carrier phase in order to define a power spectrum of $G(t)$ as defined in (1). Conditions (3) and (4) mean that the real and imaginary parts of $F(t)$ have the same autocorrelation, i.e. the same power spectrum, and that their cross-correlation is an odd function, i.e. their cross-power spectrum is purely imaginary. This last condition also implies $R_{AB}(0)=0$, i.e the real and imaginary part are uncorrelated when sampled at the same time. It turns out that for many practical problems these conditions are at least approximately satisfied.

Given that conditions (3) and (4) are satisfied, the autocorrelation function of $G(t)$ is given by

$$R_G(\tau)=\frac12\text{Re}\{R_F(\tau)e^{j\omega_0\tau}\}\tag{5}$$

and its power spectral density is given by the Fourier transform of (5):

$$S_G(\omega)=\frac14\left[S_F(\omega-\omega_0)+S^*_F(-\omega-\omega_0)\right]\tag{6}$$

where $S_F(\omega)$ is the power spectral density of $F(t)$.


It is a programming question, not a signal processing question - keep that in mind, and next time use the StackOverflow.

MATLAB rand function returns uniformly distributed variable on interval $[0,1]$.
C++ rand function returns integer between 0 and RAND_MAX, additionally it is not very 'random' function. You would also have to normalize its output. For more info please refer here: C++ rand() between 0 and 1. Although I think that you should use the uniform_int_distribution function, providing you can take the advantage of C++11 standard.

I believe that you should use OpenCV after all, instead of playing with new standards and some tricky functions. It contains the RNG class that you can easily use for generating of uniformly distributed sequence. Here is the example:

RNG rng;
float r_num = rng.uniform(0, 1);


shannon interpolation in image processing

The error source is the kernel function, which should be a multiplication of two 1D sincs (without rotational symmetry) instead of the Sombrero function, which is a sinc function with rotational symmetry.


understanding a signal from its mathematical formula

What is misleading is that you ought to find the values of $x_1[n]$ for all $n$.

For the first term, assuming convention $\cdot^0=1$, you get $\delta[n]+\delta[n-2]$. For the second one, all terms except for $k=2$ cancel each other, so you only gain another $\delta[n-2]$. Then, as you guessed, the last term (without the initial minus sign), is 1 for $n=3,4$, and zero elsewhere.


Can the order of filtering and downsampling be exchanged?

In Short: The order does matter.

Detailed:
Downsampling by a factor of 2 without filtering, will cause aliasing in frequencies of 500Hz and more.
Since your signal is bandlimited to 600Hz, the frequency range of 400Hz-500Hz will be corrupted.
So (theoretically) the frequency range 0Hz-400Hz will be the same with or without filtering.
If you were only interested in this range, it would have been possible to avoid the filtering.


How to prove relationship between mutual information and differential entropy?

Your derivation is correct. Your second integral can also be written as

$$
\int_{x,y}f(x,y)\log_2(f(y))dxdy=\int_y\log_2(f(y))\left(\int_xf(x,y)dx\right) dy
$$

where the integral in parenthesis is equal to

$$
\int_xf(x,y)dx=\int_xf(x|y)f(y)dx=f(y)\int_xf(x|y)dx
$$

You can understand this integral in two ways:

Like your idea: $\int_xf(x|y)dx=1$, which means that regardless of what value $y$ has, $x$ must have some value, so the integral over all possibilities is 1. You can understand that also by seeing $f(x|y)$ as some distribution for $x$ (which is parametrized by $y$). And, the integral over the support of a distribution needs to be $1$.
$\int_xf(x,y)dx=f(y)$ is the marginalization of $f(x,y)$ regarding the variable $x$. With the integration, you remove the influence (i.e. marginalize out) of the variable $x$ from the probability density.


Getting error in PSNR calculation using Matlab

Why this is happening?

The answer is quite simple and the error message exactly tells you what is wrong: Apparently you do not have the function psnr() available on your system.

This might have many reasons. Most probably:

You are missing a Matlab toolbox that provides this function
You are missing a Matlab script that provides this function

If you have the script, then Matlab does not know about it. Check your Matlab path settings


I assume that by "inverse Z-transform" you simply mean the impulse response of the filter.

That's much easier done in the time domain directly. The first section has an impulse response of
$$
h_1[n] =\begin{cases}
1 & \text{ if } n= 0\\
-1 & \text{ if } x= K\\
0 & otherwise
\end{cases}
$$
So basically a positive spike at n=0 and a negative spike at $n = K$

The integrator picks up the first spike and keeps it until you hit the negative spike, which cancels the positive one and everything is zero again. So the resulting impulse response is a rectangular window
$$h_2[n] =\begin{cases}
1 & \text{ if } 0 \leq n < K \\
0 & otherwise
\end{cases}$$

The second half of the filter has exactly the same structure as the first half and so it produces a rectangular window of length $L$.

The overall impulse response is then simply the convolution of both rectangular windows. This will be a triangle for $K=L$ and a trapezoid otherwise. The ramps of the trapezoid have the length of the smaller of K and L and the length of the flat part is the difference between K and L. Total length is $L+K-1$

Assuming K is the smaller one we could write this as

$$h[n] =\begin{cases}
n+1 & \text{ if } 0 \leq n < K \\
K & \text{ if } K \leq n < L \\
K-1-(n-L) & \text{ if } L \leq n < L+K-1 \\
0 & otherwise
\end{cases}$$


Wavelets can be useful feature detectors, and one of their standard advantage is to remove some trend or background effects, due to the vanishing moment property.

So first, this can act, sometimes, as a very low level non-linear effect remover.

Second, one of the successes of wavelets is the nonlinear approximation, where coefficients are selected upon their (highest) amplitude. The shrinkage/thresholding approach is nonlinear, as the shrinkage of a sum of two signals is not the same as the sum  of the shrinkage of two signals.

Third, more generally, wavelet coefficients are often used as the input of nonlinear predictors, for instance neural networks.


Is there a generalized method to get the input with the given output and the impulse response?

The output signal of an LTI system is periodic only if its input signal is periodic as well, your reasoning was right. Now, a signal periodic with a certain period $T_p$ has a line spectrum, i.e., it only contains frequency contents at $n f_0$, $n \in \mathbb{Z}$ where $f_0 = \frac{1}{T_p}$. What it basically means is that your input only excites these frequencies so that to obtain the input output relation, it only matters how the LTI system responds to these particular frequencies. You can therefore consider the impulse response of your LTI system periodic as well, it changes nothing.

That said, you can express the input-output relation entirely using Fourier series. If your input signal $x(t)$ has Fourier coefficients $C_{x}[\mu]$, your output signal $y(t)$ has Fourier coefficients $C_{y}[\mu]$ and your system is characterized by an impulse response $g(t)$ with Fourier coefficients $C_g[\mu]$ then, the input output relation is

$$y(t) = x(t) * g(t) \quad \Leftrightarrow \quad C_y[\mu] = C_g[\mu] \cdot C_x[\mu],$$

where * represents circular convolution. Therefore, if you are given $y(t)$ and $g(t)$, you can find $x(t)$ via:

Obtaining the Fourier coefficients $C_y[\mu]$ of $y(t)$
Obtaining the Fourier coefficients $C_g[\mu]$ of the (periodified) $g(t)$
Obtaining the Fourier coefficients $C_x[\mu] = C_y[\mu]/C_g[\mu]$
Finding $x(t)$ via its Fourier series expansion using $C_x[\mu]$.

Let me elaborate shortly on the difference between Fourier transforms and Fourier series: The Fourier transform takes a continuous time signal $g(t)$ and returns its amplitude density $G(\omega)$ for every $\omega$. You find it via $$G(\omega) = \int_{-\infty}^\infty g(t) {\rm e}^{-\jmath \omega t} {\rm d}t.$$ On the other hand, the Fourier series is defined for periodic functions only. It takes a periodic time signal and returns coefficients belonging to the integer multiples of its base frequency. You find the Fourier coefficients via $$C_g[\mu] = \frac{1}{t_p} \int_{t_p} g(t) {\rm e}^{-\jmath 2\pi \mu t/t_p} {\rm d}t.$$ if $g(t)$ is $t_p$ periodic. These two definitions may look very similar. Besides normalization it may seem that $C_g[\mu]$ is the same as $G(\omega)$ at $\omega = \mu \cdot \frac{2\pi}{t_p}$. However, this is only true if $g(t)$ is either actually periodic or at least time-limited to the interval $t_p$ you are considering (which is true in your case!) since in this case, periodifying it does not change it at all. If it is not limited to $t_p$, the Fourier series integral limited to the interval $t_p$ would truncate the function which would lead to wrong results. Therefore, one would have to periodify $g(t)$ first in this case to get correct results.


Of course you have to obey Shannon-Nyquist when sampling anything.

But, of course, sound cards were designed by engineers, so they have analogue filters that avoid aliasing at the sampling rates the cards offer.

So, if you're using a sound card, any sampling rate it supports should be largely alias-free. Use the sampling rate sufficient for your signal of interest and comfortable for your application and available from the sound card hardware. In effect, that means you can chose between 44.1, 48 and 96 kHz on most sound cards.


Applying a Filter Kernel Defined by a Continuous Equation in the Frequency Domain

Given an Image $ I \in \mathbb{R}^{m \times n} $ I would solve it as following (Using the Convolution Theorem):

I would sample the function in Fourier Domain into a grid of $ m \times n $.
I would apply DFT transform on the image to get the $ m \times n $ representation in the Fourier Domain.
I would apply Element Wise multiplication between the 2 $ m \times n $ arrays.
I would apply Inverse DFT on the result to get the filtered image in the spatial domain.

Pay attention that this will result in Cyclic Convolution and not Linear Convolution (Which requires more padding).


Does sampling in the frequency domain cause time-domain aliasing?

The answer is: yes, sampling in the frequency domain causes aliasing in the time domain, exactly like the dual case: sampling in the time domain causes aliasing in the frequency domain.

There are many ways to see this. One standard way is to sample the discrete-time Fourier transform (DTFT) of a discrete-time signal by multiplying it with a Dirac comb and see what this does to the time domain signal. Here, however, I would like to show it in another way, namely by taking samples of the DTFT and asking which periodic signal has these samples as its discrete Fourier transform (DFT).

Let $Y(e^{j\omega})$ be the DTFT of the sequence $y[n]$ (note that $y[n]$ is not necessarily of finite length):

$$Y(e^{j\omega})=\sum_{n=-\infty}^{\infty}y[n]e^{-jn\omega}\tag{1}$$

Now we take samples of $Y(e^{j\omega})$ at $\omega_k=2\pi k/N$, $k=0,1,\ldots,N-1$, and derive the time domain sequence with a DFT equal to these samples:

$$\begin{align}Y_k\doteq Y(e^{j2\pi k/N})&=\sum_{n=-\infty}^{\infty}y[n]e^{-j2\pi nk/N}\\&=\sum_{l=-\infty}^{\infty}\sum_{n=0}^{N-1}y[n+lN]e^{-j2\pi nk/N}\\&=
\sum_{n=0}^{N-1}\sum_{l=-\infty}^{\infty}y[n+lN]e^{-j2\pi nk/N}\\&=\sum_{n=0}^{N-1}\tilde{y}[n]e^{-j2\pi nk/N}\tag{2}\end{align}$$

where $\tilde{y}[n]$ is an aliased version of the original sequence $y[n]$:

$$\tilde{y}[n]=\sum_{l=-\infty}^{\infty}y[n+lN]\tag{3}$$

Note that the right-hand side of $(2)$ is the DFT of $\tilde{y}[n]$.
Eq. $(2)$ shows that the samples $Y_k$ of the DTFT $Y(e^{j\omega})$ are the DFT coefficients of $\tilde{y}[n]$, which is an aliased version of $y[n]$. Consequently, sampling in the frequency domain causes time-domain aliasing, which is obviously the dual of the well-known aliasing in the frequency domain caused by sampling in the time domain.


Calculating filter parameters from system response

Just math I think. Let's substitute g = ALPHA -1, which makes the algebra a little easier. Then you need to solve
$$\left | \frac{1+g}{1+g\cdot e^{-j\omega }} \right |^{2}=0.5$$
That means
$$\frac{1+g}{1+g\cdot e^{-j\omega }}\cdot \frac{1+g}{1+g\cdot e^{+j\omega }} = 0.5$$
Multiplying this out yields
$$\frac{(1+g)^{2}}{1+2\cdot g\cdot cos(\omega) + g^{2} }= 0.5$$
We can put this into standard quadratic equation form:

$$g^{2}+(4-2\cdot cos(\omega))\cdot g + 1 = 0$$

Solve, plug in the numbers and you get
$$g_{1} = -1.0744, g_{2} = -0.93074$$

The first solution doesn't make sense since it would be a pole outside the unit circle, so the second is correct and we get the expected

$$\alpha = 0.069262$$


Is there any difference between using convolution and correlation for finding edges with Sobel?

You're correct, convolution is the same as correlation except that 2D Sobel filter kernel is rotated for 180 degrees.

Generally speaking for spatial filtering using correlation or convolution is a matter of personal/domain preferences/convention.