How to compute a sweeping cutoff filter in python? I would like to compute a sweeping cutoff filter.
I wish to use SciPy.
I sliced my wave into x samples slices, in order to apply a different cutoff to each such created "frame".
The problem is that the filter, when applied to a slice filters the "clic" at the start of the slice, which creates an artefact for regathered result.
I searched stackexchange all over for sweeping filters questions and found nothing.
Any Idea on how to get a smooth sweeping without clic?
    #for filter
    from scipy.signal import butter, buttord, lfilter, freqz
    from scipy.signal import freqs

    def truncate_if_outside(data, bitwidth):
        if bitwidth==2:
            max_sample_value=32767
            min_sample_value=-32767
        else:
            print "Only 16 bit accepted by truncate_if_outside function"

        for i in range(len(data)):
            if data[i]>max_sample_value:
                data[i]=max_sample_value;
            if data[i]<min_sample_value:
                data[i]=min_sample_value;
            if math.isnan(data[i]):
                data[i]=0
            data[i]=int(data[i])
            if math.isnan(data[i]):
                data[i]=0
        return data

    def lowpass_filter(integer_data, sample_framerate, sample_bitwidth,threshold):#works with cutoff from 400 and higher

        def butter_lowpass(cutOff, fs, order):

            nyq = 0.5 * fs
            normalCutoff = cutOff / nyq
            b, a = butter(order, normalCutoff, btype='low', analog = False)
            return b, a

        def butter_lowpass_filter(data, cutOff, fs, order):
            b, a = butter_lowpass(cutOff, fs, order=order)
            y = lfilter(b, a, data)
            return y

        result=butter_lowpass_filter(integer_data, threshold, sample_framerate, order=4)#what is order?
        result=truncate_if_outside(result, sample_bitwidth)
        return result

    def sweep_filter(integer_data, start_freq, end_freq, sample_framerate, sample_bitwidth, slice_width):
        num_slice= int(len(integer_data)/ float(slice_width))
        print "nombre de tranche dans cette durée :",num_slice
        new_data=[]
        diff_cutoff=end_freq-start_freq
        for i in range(num_slice):
            current_cutoff=start_freq+((i/float(num_slice))*diff_cutoff)
            print "cutoff",current_cutoff
            print "start:",i*slice_width
            offset_first_sample=(i*slice_width)-1
            offset_last_sample=((i*slice_width)+slice_width)-1

            slice_data=integer_data[offset_first_sample:offset_last_sample]
            slice_data=lowpass_filter(slice_data, sample_framerate, sample_bitwidth, current_cutoff)
            for j in range(len(slice_data)):
                new_data.append(slice_data[j])
        return new_data

The here filtered signal is a "smoothed noise" which means a random numbers smoothed with exponential curves. There should be no gap in it.

To create a smooth sweeping filter or envelope filter, you can process as in video by "frame". You plan a 30% overlap of each frame with a crossfade.
You need to generate two tracks and mix them together to do that.


How to calculate the plane at infinity? either in Matlab or Python I tried to implement this function to compute the plane at infinity $n_\Pi$ from the equation:

where $w^{-1}$ is the image of the absolute conic IAC, $M_j$ is the Projection Matrix of image $j$,
but I'm not sure if this implementation is right, because I don't have a good idea about symbolic expressions and how they are solved.
so I just mimicked other codes, but I'm still confused, because when I used this part to calculate the intrinsic matrix, the results of the intrinsic matrix weren't converging.
can you please tell me if there are issues in my implementation, and how to fix it please?

The whole code, and data is here

Thanks in Advance

function P_infty = PlaneAtInfinity(K,P)

% Symbolic variables

X =  sym('X', 'real');
Y =  sym('Y', 'real');
Z =  sym('Z', 'real');
L2 = sym('L2','real');

N = [X; Y; Z];
W_inv = K*K';
% Quadric
Q = [W_inv       ,   (W_inv * N)    ;
    (N' * W_inv) , (N' * W_inv * N)];

% Autocalibration equation
Calibration = P * Q * P';
% func = cross(W_inv, Calibration);
% solution = Optimize(zeros(3,3),@func);
% solve linear equations
solution = solve(Calibration(1, 1) == (L2 * W_inv(1, 1)), ...
    Calibration(2, 2) == (L2 * W_inv(2, 2)), ...
    Calibration(3, 3) == (L2 * W_inv(3, 3)), ...
    Calibration(1, 3) == (L2 * W_inv(1, 3)));

P_infty = [double(solution.X(1));...
    double(solution.Y(1));...
    double(solution.Z(1))];
end

and the main script:

P = [4.35310807409356, 1.44651724424134, 36.0849770394682, 6.6274893804074;
22.9809618064503, 8.56379922408173, -2.94837061198918, 36.6987313462383;
-0.00438517470946907, -0.00160149495653558, 0.241477646472821, -0.00696180145320084];

K = [870.258271610840, 0, 279.100746891319;
0, 812.178628533806, 260.930050222247;
0, 0, 1];
PlaneAtInfinity(K,P)

The plane at infinity has to be calculated as following:
function plane = computePlaneAtInfinity(P, K)

%Input
%   P - Projection matrices
%   K - Approximate Values of Intrinsics
%
%Output
%   plane - coordinate of plane at infinity

% Compute the DIAC W^-1
W_invert = K * K';

% Construct Symbolic Variables to Solve for Plane at Infinity
% X,Y,Z is the coordinate of plane at infinity
% XX is the scale
X = sym('X', 'real');
Y = sym('Y', 'real');
Z = sym('Z', 'real');
XX = sym('XX', 'real');

% Define Normal to Plane at Infinity
N = [X; Y; Z];

% Equation of Dual Absolute Quadric (DAQ)
Q = [W_invert, (W_invert * N); (N' * W_invert), (N' * W_invert * N)];

% Select Any One Projection Matrix
M = P(:, :, 2);

% Left hand side of the equation
LHS = M * Q * M';

% Solve for [X, Y, Z] considering the System of Linear Equations
% We need 4 equations for 4 variables X,Y,Z,XX
S = solve(LHS(1, 1) == (XX * W_invert(1, 1)), ...
          LHS(1, 2) == (XX * W_invert(1, 2)), ...
          LHS(1, 3) == (XX * W_invert(1, 3)), ...
          LHS(2, 2) == (XX * W_invert(2, 2)));

plane = [double(S.X(1)); double(S.Y(1)); double(S.Z(1))];

end
```


best meathod to lock onto a cdma signal? I've been looking at 2 different methods to acquire a known prn code in a CDMA signal. The first trys Xor prn code chip to I sign data. Then average the result over a time period to determine the data bit. I have not had reliable success with this method. The second which I am going to try performs a real multiply of a the code chip and the sign then adds or subtracts the I data magnitude value over a time span.
Is one method known to be more reliable than the other? I'm typing on phone so I can't post any coding right now. If it would help I can update this question later.
Thanks

CDMA signals typically have a pilot that is sent along with the data.  This pilot is a known chip sequence that is transmitter over and over again.  Receivers can use the pilot to lock onto the signal by correlating the received signal with the expected pilot.

The autocorrelation is typically done by FFT'ing the received signal, multiplying it by the FFT'ed expected pilot, and IFFT'ing the result.  This is done both for computational efficiency and because the circular correlation that you get with FFT's is actually helpful in this case because the pilot is repeated.  Thus you should be able to "catch" the pilot no matter which phase it starts at in your captured signal.

Once you have locked onto the pilot you know the signal's phase and it should be pretty straightforward at that point to despread the signal and recover the data, providing you know the data's spreading sequence.


Microscopy crystal particle object length detection Hi I  try to detect the particle length size from this image

and I'd like to get the length size of the particle Like this below. and get the histogram of the size of the particle using centroid and detect the line like this below:

I have try to using watershed. But i still low knowledge on this code. What i want to ask are:
1. How to get segmentation like this picture as an example below so I can evaluate my segmentation.. I dont know how to display it (I will share my code)

How to separate the particle that unite with other particle. I got this contour

Here it is my code:

from __future__ import print_function
from imutils import perspective
from imutils import contours
import numpy as np
import argparse
import imutils
import cv2
import scalar
import sys
from scipy.spatial import distance as dist
import scipy.ndimage as ndi
import matplotlib.pyplot as plt
from skimage import filters
from sklearn import cluster
from scipy import ndimage
from skimage.feature import peak_local_max
from skimage.morphology import watershed
from scipy.ndimage import label

img = cv2.imread("C:\\Users\\adiyu\Pictures\\snipping\\c3smeasure.png",
cv2.IMREAD_COLOR)
kernel = np.array([[-1,-1,-1],[-1,9,-1], [-1,-1,-1]])
kernel2 = np.ones((3,3),np.uint8)

img=cv2.filter2D(img,-255, kernel)
img = cv2.medianBlur(img,3)
img_np = np.array(img)
img_np_rgb = cv2.cvtColor(img_np,cv2.COLOR_RGBA2RGB)

#Kmeans cluster using Opencv
Z = img.reshape((-1,3))

# convert to np.float32
Z = np.float32(Z)

# define criteria, number of clusters(K) and apply kmeans()
criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 10, 1.0)
K = 3
ret,label,center=cv2.kmeans(Z,K,None,criteria,10,cv2.KMEANS_PP_CENTERS)
# Now convert back into uint8, and make original image
center = np.uint8(center)

res = center[label.flatten()]
res2 = res.reshape((img.shape))
gray= cv2.cvtColor(res2,cv2.COLOR_BGR2GRAY)
gray = cv2.fastNlMeansDenoising(gray, None, 6,7,21)

cv2.imshow('res2',res2)
cv2.imshow('gray',gray)
thresh = cv2.threshold(gray, 100, 255, cv2.THRESH_BINARY_INV +
cv2.THRESH_OTSU + cv2.THRESH_OTSU)[1]
cv2.imshow('thresh',thresh)
opening = cv2.morphologyEx(thresh, cv2.MORPH_CLOSE, kernel, iterations=2)
cv2.imshow("opening", opening)
closing = cv2.morphologyEx(opening, cv2.MORPH_CLOSE,kernel2, iterations=2)
cv2.imshow("closing", closing)
fg = cv2.erode(closing,None,iterations = 2)
cv2.imshow("fg", fg)

# sure background area
sure_bg = cv2.dilate(opening,kernel,iterations=3)
cv2.imshow("sure_bg", sure_bg)
# Finding sure foreground area
dist_transform = cv2.distanceTransform(opening,cv2.DIST_L2,5)
ret, sure_fg = cv2.threshold(dist_transform,0.1*dist_transform.max(),255,0)

# Finding unknown region
sure_fg = np.uint8(sure_fg)
unknown = cv2.subtract(sure_bg,sure_fg)
cv2.imshow("unknown", unknown)

# Marker labelling
ret, markers = cv2.connectedComponents(sure_fg)

# Add one to all labels so that sure background is not 0, but 1
markers = markers+1

# Now, mark the region of unknown with zero
markers[unknown==255] = 0
markers = cv2.watershed(img,markers)
img[markers == -1] = [255,0,0]

cv2.imshow("img", img)

cv2.waitKey(0)
cv2.destroyAllWindows()enter preformatted text here

any suggestions are welcome..Thank U

Does it have to use python? For exploring your workflow and pipeline for a kind of image analysis, I would suggest getting ImageJ / FiJi. Use it together with the MorphoLibJ plugin to get more options on morphological segmentation.
Once you determined your workflow you can code it in python to batch analyze.

Anyways, your workflow should look like this, regardless of which software you use:

Segment particles either by rgb-thresholding or by converting to grayscale beforehand. An adaptive method is preferred. You can also use the contours you got. The aim is to produce a binary (black/white) image, where your particles are white and everything else is black (or the other way around).
On this image you perform "distance transform watershed". Which is also doable with python.
Perform optical checks on the resulting image to confirm that particles aren't over- or undersegmented. Refine parameters like the dynamics, the distance metric, seeding, etc. if needed.
If correctly segmented perform particle analysis on the image. In ImageJ you can go "Analyze -> Set Measurements" and select variables to measure - for example Feret Diameter - Length. In python i am not experienced enough but I would go with this tutorial


Are radar resolution and detection capabilities not very tightly related? Radar detection capacity is the capacity to determine the presence or absence of targets whereas the resolution is the capacity to identify two closely separated targets.  If you can detect the target separately , you have resolved the targets.   If you resolved two targets, you have detected them too.  Thus are the detection and resolution tightly bonded if not the same?  I think, I am wrong some where and I want to know it.

Resolution is (usually) referred to as the ability to distinguish two closely separated returns, not targets.

The distinction is how far along the processing chain of the radar you are, resolution, given by the inverse of the bandwidth is often limited by the hardware (waveform generation, filters, amplifiers, sampling..) and the environment. While Radar detection is done as a further processing to distinguish between what you want to see (targets) and clutter (birds, ground, sea, trees, ..) edit: and is usually limited by the signal to noise/clutter ratio. The detection step will usually contain multiple looks.

You can therefore have say a really good resolution, but since all you see is "noise/clutter", you can't really detect the target you are looking for. Edit: Conversely you may have a really good detection algorithm that uses say excessive averaging and multiple 'looks', but if your resolution is poor, two closely spaced targets may appear as a single target.

Edit: Note that I use noise and clutter interchangeably, as it simply refers to what you don't want to see.

Edit: To simplify things: resolution will depend on the bandwidth, while detection will depend on noise level and processing.

I have attempted to sketch a generic radar system below, the resolution is given by the receiver output, while the detection algorithm will typically use multiple looks (for averaging and doppler discrimination) and use a threshold to alert the user of one or more target.


Q matrix and updating times in a Kalman filter The context of the problem is that I have several robots located remotely which give their position (x,y coordinates) every x seconds and send it to a centralized remote server. The value of the seconds, x, depends on the update frequency of the robot, so it can be 60, 70, 80 or whatever seconds.
In this server, I take these measurements (x, y coordinates of each robot) and generate a Kalman filter to predict next steps, since the measurements are given each "too much" time and I need these predicted positions.
The Kalman filter has the equations:

x[k] = x[k-1] + v[k-1]*dt
v[k] = v[k-1]

The robots are linked by pairs, and the filter of each robot (inside a pair) is updated when a measurement arrives from any of the two robots.
Here comes the problem. The frequency of update in the filters depends on the time of the robots' measurements. For instance:

Two robots giving measurements every 10s and every 25s. The filters will be initialized at t=0s and then updated at t=10s, t=20, t=25s, t=30s, etc.

Two robots giving measurements every 3s and every 5s. The filters will be initialized at t=0s and then updated at t=3s, t=5, t=6s, t=9s, t=10s, etc.

Because of this different updating times, the filters of some pairs of robots work very well, but others do not.
The matrixes Q and R stays constant at any moment. For the Q matrix, the values where given acording to the error produced by the maximum aceleration that can be reached (the problem should be here).
How can I implement a filter (a class) that works for all of the pairs independently of the updating filter times and updated measurement times? Which is the dependence between the update time of the filter and the time when a measurement arrives? I mean, how does it affect quantitatively?
Should I modify the matrixes Q, R and/or P according to the time diference?
It is an open question since I am a bit lost. Any help would be highly appreciated.

Kalman filter works in a predictor-corrector or predict-update sequence. When no new measurement is available the only thing that you can do is to predict(NOT update) the state with your assumed dynamic model. The prediction is done not only on the state vector itself but also on its covariance matrix ${\bf{P}}$ (Not the measurement covariance matrix $\bf{R}$). A popular model commonly used for model regularization (AKA model noise matrix, or the matrix $\bf{Q}$) is the Singer model from 1990 (see [1]).

In your case the transition matrix is given by $\Phi_k =  \begin{bmatrix}1 & dt \\0 & 1 \end{bmatrix}$ , so the state covariance matrix is predicted using
${\bf{P}}_{k|k-1}=\Phi_k {\bf{P}}_{k-1|k-1} \Phi_k^T + {\bf{Q}}_k$. That way the state matrix ${\bf{P}}$ "inflates" to account for the uncertainty.

[1] Singer, R. A., “Estimating Optimal Tracking Filter Performance for Manned Maneuvering Targets,” IEEE Transactions on Aerospace and Electronic Systems, vol. AES-6, no. 4, .pp. 473-483, July 1990.


A 23 kHz sine wave sampled at 48 kHz, then FFT A 23 kHz sine wave sampled at 48 kHz with 100 samples. Then do a 64 point FFT of the sampled signal. Why is the FFT of the signal not a delta function at 23 kHz. Or why does the  FFT also have values at other frequency?

I think it has something to do with 23kHz / 48 kHz * 64 samples not being an integer. I just couldn’t think the logic through. Please help me! Thank you!!

Edited: sorry guys, i meant 23 kHz

Since you used 2.3k twice, I am assuming you meant that with the original 23k value.

Converting your frequency:

$$ 2.3k \frac{cycles}{second} \cdot \frac{ 64 \frac{samples}{frame} }{ 48k \frac{samples}{second}  }  \approx 3.07 \frac{cycles}{frame}$$

Because this is not an integer value, the peak DFT value will be at bin 3 and there will be values in all the other bins, which is commonly called "leakage".  Being a real valued signal, the upper half of the DFT will be a mirror image with complex congjugate values, so there will also be a peak at bin 61 (64-3).

Your 100 sample figure is irrelevant if you are only using 64 of them for the DFT.

Signals have to be a whole number of cycles per frame in order for the DFT not to have leakage.  The other answers presume you understand that the DFT is a dirac delta, windowed, convolved version of the continuous FT.  This is not an easy way to understand it, but it is how it is usually taught.  It also leads to the misconception that the leakage is described by what is known as the sinc function $ \frac{\sin(d)}{d} $.  That is not true.  The fall off for a complex valued signal is described by what is known as the aliased sinc function $ \frac{\sin(Nd)}{N\sin(d)} $.  Another name is "Dirichlet function".  A real valued signal is actually the sum of two complex signals, so that function doesn't apply either.

You can find the exact formula for the bin values, the peak and "leakage", in my blog article "DFT Bin Value Formulas for Pure Real Tones".

You can also find a much more intuitive explanation of why leakage occurs in another one of my blog articles "DFT Graphical Interpretation: Centroids of Weighted Roots of Unity".

Hope this helps,

Ced


What is a 'time-variant' signal? I have found in the literature a "time-variant signal" term , but I could not understand what that's mean. I know that there is "time-variant systems", but this is the first time, which I see a time-variant signal term. Any one has an idea what this term means ?   Some suggests that this term is related to non-stationary signal, which it's spectrum changes over time. Is this right?

P.S: I have added the  EEG tag because this term, is used widely to describe the EEG signals

As I've mentioned in a comment, time variant signal is a term which is not commonly used in signal processing, and I also couldn't find a clear definition. After browsing a few texts which use the term (usually written by people not specialized in signal processing), I've come to the conclusion that by time variant signals they actually mean non-stationary signals. The latter being a clearly defined term (check the wikipedia page on stationarity). In short, it means that the signal's properties, such as its mean and variance, change over time. And this is indeed the case with EEG signals, just as it is the case with speech signals.


Fractional powers of complex numbers (DSPrelated computation) I am puzzled by this computation: https://www.dsprelated.com/showarticle/754.php (c.f. quote)

  Raising $ i $ to integer powers results in traversing the unitcircle in the same number of quarter turns. The next question that arises naturally is if the pattern also applies to fractional values as well. Consider $ \sqrt{i} $. Suppose it is $ a + b \cdot i $.
  $$ ( a + b \cdot i )^{2} = i \\ a^{2} + 2ab \cdot i + b^{2} \cdot i^{2} = i $$
  $$ ( a^{2} - b^{2} ) + ( 2ab ) \cdot i = 0 + 1 \cdot i $$ In order for two complex numbers to be equal, both the real parts and the imaginary parts must be equal. $$ a^{2} - b^{2} = 0 \text{ and } 2ab = 1 $$
  $$ a = ± b \text{ and } a = { 1 \over {2b} } $$
  Solving for a,b being real yields two solutions: $$ (a,b) = \left( { \sqrt{2} \over 2 }, { \sqrt{2} \over 2 } \right) \text{ or } \left( - { \sqrt{2} \over 2 }, - { \sqrt{2} \over 2 } \right) $$

For me, it does not make sense to say that $a=\pm b$ AND $2ab=1$, this is just a contradiction... either $a=\pm b$ OR $2ab=1$ therefore I would say that $(a+bi)^2 \neq i$ .... How can they say that it is? What is the logic behind that...?

That's my article, I'm glad you are reading it.

A fundamental rule for two complex numbers (corresponding to a point on a plane) is that they are only equal if both the real and imaginary parts (corresponding to x,y coordinates) are equal.  Therefore the complex equation

$$ (a^2 - b^2) + 2abi = i = 0 + i $$

is equivalent to saying the point $(a^2 - b^2,2ab)$ is the same point as $(0, 1)$

The real parts must be equal AND the imaginary parts must be equal.

Graphically, you have the intersection of a hyperbola intersecting crossed lines going through the origin.  Substitute $x$ and $y$ for $a$ and $b$ if that helps you see it.

The important concept is that an intersection point is halfway along the unit circle, and the other one is directly opposite (the negative of it.)  This is true for taking the square root of any number on the unit circle.  One square root will be halfway and the other the negative of it.

And that is part of the exponential nature of the unit circle.  Radians is just one possible scale.  Quarter turns is another, whole turns another.

Expanding on this, the solutions to $1^{1/N}$ are $N$ evenly spaced points around the circle.  They are called the Roots of Unity.  You may know the first root better as $e^{i2\pi/N}$ which is right there in the DFT definition:

$$ X[k] = \sum_{n=0}^{N-1} x[n] \left(e^{i2\pi/N}\right)^{-nk} $$

$$ a^2 - b^2 = 0 $$

$$ a^2 = b^2 $$

$$ a = \pm b $$

And

$$ 2ab = 1 $$

For $a=b$

$$ 2a^2 = 1 $$

$$ a = \pm \frac{\sqrt{2}}{2} = b $$

For $a=-b$

$$ -2a^2 = 1 $$

No real solutions.


Extracting color information from lit regions to adjust shadow regions, HSV space If I wanted to match the color of a shadow region to the color it is supposed to be in the lit region, would this, purely color information, be preserved in the hue of the HSV space of the image?

Probably not – the observed color of a material is a result of the composition of light that illuminates the spot, and what of that the material reflects.
For example, if you have a bright light with strong red/yellow component ("warm light"), and you have ambient/scattered light at a lower level that's more or less white, then a green spot in the bright area will have a different tint than in the shadow.
(note: a shadow isn't perfect darkness, otherwise everything there would be the blackest of blacks.)
So, you need have a model of at least your main light and how ambient light, eg. from scattering that on a not-quite-perfectly-white wall, changes the spectrum of that, and of the actual color of the material (which you can't see in either, unless one of these is a reference white light source...).


How does a Hermitian FFT work in Numpy? Say, I create a Hermitian complex signal using,

import numpy as np
t = np.arange(-4, 4)
z = np.exp(1j * t)

Here z should be a complex signal with Hermitian Symmetry, as you can see below.

In [2]: t
Out[2]: array([-4, -3, -2, -1,  0,  1,  2,  3])

In [3]: z
Out[3]:
array([-0.65364362+0.7568025j , -0.98999250-0.14112001j,
       -0.41614684-0.90929743j,  0.54030231-0.84147098j,
        1.00000000+0.j        ,  0.54030231+0.84147098j,
       -0.41614684+0.90929743j, -0.98999250+0.14112001j])

But when I take the Fourier Transform of this signal, I don't get a real spectrum,

In [8]: np.fft.fft(z)
Out[8]:
array([-1.38531768+0.7568025j, -7.02599565+0.7568025j,
        2.57935201+0.7568025j,  0.93952139+0.7568025j,
        0.41344309+0.7568025j,  0.08151870+0.7568025j,
       -0.22205190+0.7568025j, -0.60961892+0.7568025j])

However, when I take the Hermitian FFT using the hfft function, I get

In [10]: np.fft.hfft(z[:5])
Out[10]:
array([-1.38531768, -7.02599565,  2.57935201,  0.93952139,  0.41344309,
        0.0815187 , -0.2220519 , -0.60961892])

This is the real component of the result I got using the regular fft function.
I don't understand what I am doing wrong here. Shouldn't fft be giving a result with a zero or almost-zero imaginary component when the input is Hermitian?

I feel like I am doing something wrong here regarding how the signal is presented to fft.

You are facing the practical problem of representing a noncausal signal in a programs such as Matlab, Octave, or Python. These programs, in their default, assume that all represented sequences begin at time-index $n=0$.
Therefore, a Hermitian symmetric sequence ( $x[n] = x^*[-n]$ for all $n$ ) ,  which begins at $n=-4$, and ends at $n=4$, will lose this symmetry when interpreted in those programs. The effect is a shift of the signal to the right by $4$ for your particular example.
Hence, what you have compouted is not the FFT of the original symmetric sequence $x[n] = e^{j n}$, but a shifted, non-symmetric, sequence $w[n] = x[n-4]$. Using the DTFT / DFT time-shift property, we see what is actually computed as :
$$ W(e^{j\omega}) = e^{-j4\omega} X(e^{j\omega}) \tag{1} $$ and
$$ W[k] = e^{-j4 \frac{2\pi}{N} k} X[k] ~~~~~ , ~~\text{ for } ~~~k=0,1,...,N-1 \tag{2}$$
Where $N$ is the length of the sequences $x[n]$ and its DFT $X[k]$.
Therefore, eventhough the DTFT / DFT of the desired sequence $x[n]$ has zero imaginary part, as an Hermitian symmetric signal would have, the practically computed $W[k]$ has nonzero imaginary part.
If you wish, you can obtain the non-zero phase samples of the desired DFT $X[k]$ of $x[n]$, from the computed DFT $W[k]$ of $w[n]$, with the following:
$$ X[k] = e^{j 4 \frac{2\pi}{N} k} W[k] ~~~~ ,~~\text{ for } k=0,1,...,N-1 \tag{3}$$
The following MATLAB code demonstrates the issue:
% Demonstrates the Hermitian symmetry and its handling in Matlab/Octave/Python
% environments...

a = -4;
b = +4;
n = [a:b];
x = exp(1j*n);      % Samples of Hermitian symmetrix input: x[k] = x[-k]*

figure,stem( imag( exp(-1j*2*pi*a*[0:(b-a)]/(b-a+1)).*fft(x) ));

The output is:


FastICA - Purpose of centering and whitening the Data I'm learning about ICA, but dont know whats the purpose of centering the data (making each component have zero mean)and whitening the data (not sure what is this for, is it the same as PCA?), I'm following FastICA algorithm for wikipedia
FastICA

The goal of FastICA is to rotate your data (unitary transform) so that each axis looks as non-Gaussian as possible.  Gaussian data still looks Gaussian when you rotate it.

If you don't "sphere" the data, all the algorithm can really do is rotate the whole block to one axis.  By bringing the mean to zero (centering), and normalizing the variance in all directions (whitening), you give the algorithm freedom to rotate in all directions.

This lets FastICA find the rotations that correspond to non-gaussian data.  All of the rotations it tries will leave the Gaussian components of your data as zero-mean and unit-variance, but non-gaussian components can be optimized.


What are the eigenvalues of the 8 point DFT matrix? I  know that the eigenvalues for 4 point DFT matrix can be found from $F_4^4=I$. Is this also valid for 8, 16 and higher orders? For example with  8 points, will it be $F_8^8=I$ ? If not, how  can I can compute them?

The eigenvalues belong to the same set of quartic roots of unity verifying $\lambda^4=1$, whatever the order of the DFT.

For more details on their multiplicity, you can read:
Eigenvectors and Functions of the Discrete Fourier Transform, 1982, Dickinson and Steiglitz (online).


Performance bounds for low-pass filters Precisely how does the performance of an FIR filter depend on its design criteria?
To illustrate, here is a specific question that probably does not have a convenient answer:

The Parks-McClellan algorithm converges to a FIR filter optimally aligned to design parameters. With it one can design a low-pass filter with minimum pass- and stop-band ripple:

Take $\mathcal{S}_n$ to be the collection of symmetric FIR filters of order $n$:
\begin{align}
\mathcal{S}_n := \left\{(h_1,\dots,h_n)\in \mathbb{R}^n: h_k=h_{n-k+1} \text{ for each }k\right\}.
\end{align}
Write the frequency response of $h\in \mathcal{S}_n$ as $\hat{h}(\omega)=\sum_k h_ke^{-j2\pi \omega k}.$
Specify a "stopband weight" $\alpha > 0$, cutoff $\omega>0$ and transition width $\varepsilon>0$. Define filter mismatch as follows:
\begin{align}
    \gamma_n(h; \alpha,\omega,\epsilon)\!
:= \min_{f\in [0,\omega]} \left|1-|\hat{h}(f)|\right| + \alpha\cdot \max_{f\in[\omega+\epsilon, 2\pi]} |\hat{h}(f)|.
\end{align}
The Parks-McClellan algorithm with parameters $(\alpha, \omega,\varepsilon)$ converges to an $h_{best}$ which minimizes $\gamma_n(h; \alpha, \omega, \varepsilon)$ over all $\mathcal{S}^n$.

What estimating expressions depending on $(n, \alpha,\omega, \varepsilon)$ provide tight upper and lower bounds for $\min_h \gamma_n(h;\alpha, \omega,\varepsilon)$ ?
Alternatively, fixing $n,\ (\alpha,\omega, \varepsilon)$  and $\Delta > 0$, are there bounds for $\min_h \max_{f\in[\omega+\epsilon, 2\pi]} |\hat{h}(f)|$ among those $h\in \mathcal{S}_n$ that have $\min_{f\in [0,\omega]} \left|1-|\hat{h}(f)|\right| \leq \delta$

As far as I know there are no known tight upper and lower bounds, but there's a quite accurate empirical formula for estimating the required filter order $M$ given maximum passband and stopband deviations $\delta_1$ and $\delta_2$, and the width of the transition band $\Delta\omega$:
$$M=\frac{-10\log_{10}(\delta_1\delta_2)-13}{2.324\;\Delta\omega}\tag{1}$$
Formula $(1)$ can be found as Eq. $(7.117)$ in Oppenheim and Schafer's textbook Discrete-Time Signal Processing (3rd ed, p.570). According to the authors, this formula is due to Kaiser, who developed it based on experimental results published by Herrmann, Rabiner, and Chan.
Also Matlab's firpmord.m uses a similar but more complex formula to estimate the required filter order.


Correct comparison of FFT noise with FFT signal I have a sensor which I am using to track a strong but noisy varying signal. When the signal is present there is a lot of energy in the FFT bins, when there is no signal the is much less energy - but still not zero.

I am wondering what is the correct way to use the no-signal FFT data (at startup I know there is no signal) to determine whether there is actually a signal at a later time.

Currently I am creating a noise reference by taking an FFT of the background noise and averaging it across all the bins. I am then using that reference via simple comparison to determine whether a bin contains a signal or not.

This works ok, but not perfectly and I am wondering whether calculating the Power Spectral Density of the noise and then comparing that bin-by-bin with the PSD of the signal would give better results. I would plan to scale the PSD as per https://holometer.fnal.gov/GH_FFT.pdf to eliminate the effects of the Hanning window that I am using.

I can't use ensemble averages as time resolution is important, although I could generate an ensemble average of the noise at startup if that would help.

Update: I found this paper: https://pdfs.semanticscholar.org/e59e/303e8f8fca708fe557bc36babb254ffa07f8.pdf which seems almost exactly what I want but too complex. It seems to suggest finding a noise floor by taking the minimum PSD from the spectra and then doing a direct comparison with the PSD of the signal - but clearly its got to be more compicated than that. Could I use this to calculate SNR of the peak bin and configure the minimum acceptable value?

It's interesting that there are no slam-dunk answers here - I think that's good, it means that there is no right answer and what works is going to hinge on the specifics of my implementation. Here is what I did:

I switched to measuring the power spectrum normalized for the Hanning window. This seems to give me a more stable center frequency and reduces the noise when there is no signal
I measured the noise on a per-bin basis and recorded the power spectrum and then averaged the output over a series of periodigrams. Since this reduces the gaussian noise I then multiplied by sqrt(N) where N was the number of periodigrams that I collected. I realise that this will amplify the non-random noise, but it also gives me a reasonable noise floor for the gaussian noise
I compared the noise power with the signal noise, converted to dB and allowed this to be configured by the user. 25dB works well in my case. I special-cased the situation where the noise was < 1 and made this the floor

This seems to work pretty well.


How to detect if a signal raises sharply or slowly? so how can i detect if a signal raises sharply or not ? Are there easy solutions (i suppose there are) ?

the problem is that if a signal raises sharply, i need less data to not pollute the dataset and if the signal raises slowly then i need more data. But how to determine the number of data, without knowing if it's sharp or not?

Also I need to know if it raises sharply when i'm at the top of the curve, NOT afterwards (in realtime if you want)

thanks so much

here is two examples of signals
http://khaelis.com/downloads/screenshotSharp.jpg
http://khaelis.com/downloads/screenshotSlow.jpg

Jeff

I believe that you can calculate the derivate which is the rate of signal change. In discrete signals this is obtained (approximated) by:

$$x'[k] = x[k]-x[k-1]$$

(by other words, it's the difference between consecutive samples)


Why is it necessary to implement octaves in sift I've been studying sift very hard for two weeks.
I found much materials about scale space.
It is very hard to understand scale space depthly.

What I've found and confusing things are that

'scale' in scale space means the sigma of gaussian function.
'octave' is used to represent images that are reduced.

In paper by Lowe, for each octave, images are represented. Is it not necessary? In other words, is only scale space sufficient?

I've seen many web sites as well as here, but scale space is very hard. Is there anyone who can explain in detail?

I don't know if I completely understand your question, but I will have a go at clarifying the scale space, multi-resolution ocataves and why they are important for SIFT.

To understand the scale space it is helpful to consider how you recognise images at different distances (e.g far away you may be able to distinguish the shape of a person. As that person moves closer more specific details may become apparent and thus you may recognise them). So, the spacial scale is important when attempting to recognise features in an image.

However, when we are processing an image we do not know the spacial scale at which the important features for recognition will be present. Therefore by blurring the image with a Gaussian filter, each image represents a different spacial scale. It is important to note that when a distinctive feature appears at a given scale, it will persist in images that have been smoothed with a smaller kernel. An interesting question is then if we can identify which spatial scale a given feature presents itself at. This can be done by taking the difference of two Gaussian filters. This is termed the difference of Gaussians (DOG). This is precisely what is done in the SIFT algorithm - and is thought to approximate the human visual system to some degree. SIFT achieves scale invariance by looking for extrema across different spatial bands.

It is here the multi-resolution pyramid comes into play. Because blurred images represent lower frequency information it is computationally more concise to downscale the image at each frequency band. This allows for a uniform method to explore the spacial scale of the information. Remember, the purpose of SIFT is to find features which are as invariant as possible across different scales. So being able to explore the spacial scales at which features occur is important.

I hope this helps. You may also find http://www.aishack.in/2010/05/sift-scale-invariant-feature-transform/ useful.


What type of information can be contained in I/Q (quadrature) data? I am taking in regards to, for example, the signals captured using a measuring instrument such as a spectrum analyzer. Basically would the I/Q data store all the parameters that my analyzer can analyze? What is so great about I/Q data? I would like to understand because, clearly, I do not...
Thank you.

IQ data samples can contain double the bandwidth of information within the same frequency bandwidth and/or sample rate of single channel (non-IQ or strictly-"real") data samples.

When heterodyning a real RF signal to a strictly real signal at another frequency (usually a lower IF frequency or to baseband, for easier further filtering, processing, or demodulation, etc.), the upper and lower sidebands of the original signal around the frequency of the heterodyning local oscillator will be aliased together into half the bandwidth.  The upper and lower half of an FFT of that strictly real resulting signal will contain redundant (conjugate mirrored) result bins (only half of the full complex FFT result will contain "useful"/independant information of the aliased mix).

When heterodyning a real RF signal to an IQ signal of that same bandwidth and/or sample rate, the upper sideband and lower sideband will not be aliased together, and thus can represent twice the bandwidth. An FFT of the resulting IQ data can contain independant and different data in its lower and "upper" half of the FFT's complex result bins, representing both the unmixed upper and lower sidebands of the original RF signal.  Thus twice as much frequency information in the IQ heterodyne as compared to the single channel strictly-"real" heterodyne (of the same sample rate or result channel spectral frequency bandwidth).


Why using the impulse respond of an LTI to characterize its output? We all know that every LTI, is defined by its transfer function, so we can determine the output of any system, (using common transforms). It is also known that we can determine the output of any LTI for any given input signal just by having its impulse response (which gives us an output equal to its transfer function). I don't know if I'm following the point here, is it possible to determine the impulse respond of a system without having its transfer function? If so, how is it obtained?
Can you give me a hand with this please?

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)
$$


What is a good FFT window function to reject DC? I'm using an FFT to analyze what is essentially the power envelope of a signal (see here for info on the containing project), and, since power numbers are always positive, to eliminate the DC component I'd like to use a window function that is 50/50 positive and negative, vs the usual all-positive function.

I've taken the "flat top" function, removed the a0 bias and converted it from cosines to sines, but I'm not sure that's optimal (or even meaningful).

Any suggestion?

If you're concerned with doing spectral analysis on a signal with a large DC component, and you want to suppress that DC peak, then a window function is not what you want. As some other answers noted, a highpass filter (or, viewed differently, a notch filter with the notch at zero frequency) is an appropriate solution.

To understand why, you need to think about what applying a window function does to the frequency response of each DFT output. The DFT is defined as:

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

One interpretation of how the DFT works is as a bank of filters at $N$ equally-spaced frequencies between $-\frac{f_s}{2}$ and $\frac{f_s}{2}$. Recast the sum above as follows:

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

where:

$$
x_k[n] = x[n] e^{\frac{-j2 \pi n k}{N}}
$$

So, the $k$-th DFT output is generated by first taking the input signal $x[n]$ and multiplying it by a complex exponential at frequency $\frac{-2\pi k}{N}$ to yield a downconverted signal $x_k[n]$. The resulting signal is then summed over the $N$-sample window to yield the DFT output $X[k]$. This is effectively a moving average filter (sometimes called a boxcar filter), whose impulse response can be described as:

$$
b[n] =
\begin{cases}
1,\ x = 0, 1, \ldots , N-1 \\
0,\ \text{otherwise}
\end{cases}
$$

The magnitude response of the boxcar filter can be found by taking the discrete-time Fourier transform (DTFT) of that impulse response:

$$
|H(f)| = \left|\frac{\sin\left(N \pi\frac{f}{f_s}\right)}{\sin\left(\pi\frac{f}{f_s}\right)}\right|
$$

This is a Dirichlet kernel, and is sometimes referred to as a "periodic sinc" since it looks a bit like a sinc function but repeats periodically, which a sinc does not. This expression gives the magnitude response of each DFT output, where $f$ is measured as the frequency offset from the center frequency of the respective output bin. This illustrates the spectral leakage effect; each DFT output has a frequency response that covers some continuous swath of the input signal's spectrum, not just the discrete center frequency of each output.

Now consider how things change if you apply a window function to the input signal $x[n]$ before performing the DFT:

$$
\begin{align*}
X[k] &= \sum_{n=0}^{N-1} w[n] x[n] e^{\frac{-j2 \pi n k}{N}} \\
&= \sum_{n=0}^{N-1} w[n] x_k[n]
\end{align*}
$$

With the window function in place, the downconverted $x_k[n]$ is effectively passing through an FIR filter with an impulse response described by the window function. So, the per-output magnitude response of the DFT is:

$$
|H(f)| = |W(f)|
$$

where $W(f)$ the DTFT of the window function $w[n]$. Now note that if you chose a window function that had a zero at DC and used it to premultiply $x[n]$ before the DFT, you would actually cause the unintended effect of nulling out not only DC in the resulting spectrum, but the center frequencies of every one of the DFT outputs. This is probably not what you want.

So, if you truly just want to cancel the signal's DC component, removing it via some other type of pre-processing, not time-domain windowing, is the way to go. You could use a linear highpass filter with a very low cutoff frequency or subtract the estimated mean from the signal first, for example. Choosing between these methods should be based upon what other constraints your system has.


Fourier transform of Image to identify sinusoidal sources of interference I'm a statistics grad student, and I just started getting into Digital-Image-Processing (an analogy for processing super-large contingency tables).  In the book "Digital Image Processing" by Gonzalez and Woods, I was reading chapter 2, page 94 of the third edition, and found the following image that is MAGIC to me...

They begin with a noisy image, identify some 'hot points' surrounding the center of the Fourier transform, then simply remove them, reverse the transform, and I can't help vocally demonstrate an emotional reaction while assessing that the image is now crystal clear.

Does anyone have any insight into this?  I'd like to know what are the underlying notions and methods used to filter this image.

Additionally, if any of the StackExchange users a familiar with this topic and can suggest some articles/books that I might read, I would be very grateful.

Alright. This is image processing. Sometimes its results might seem indeed magical. And quite some other times, in some books or papers, you see they discuss an algorithm or method to process a given distorted/noisy image and convert it into an unbelievably perfected form. Yet they don't provide any actual code, neither are the input and output images made available. It's therefore quite impossible to replicate what they propose to have achieved.

In this example a ridiculously simple idea is provoked. In fact it is the most basic idea of all signal procesing: Filtering! There is an original degraded image so much degraded that it might seem impossible to recover, but when it's analysed via DFT it's been seen that there're some spurious frequency components which are not to be expected from the original undistorted image. In addition it's quite interesting that those spurious frequencies are all distributed inside a circular ring of frequencies! So this is the easiest of all signal processing. Just design a circular 2D filter narrow enough to suppress all those circular band of spurious components, and leave all the rest untuoched. The result is that magic.

That being said, however, It is very unlikely that in your everyday life you would encounter such a highly structured "noise" or "distortion". It seems as if almost this distortion was just added there for the purpose of a filtering demonstration. And indeed due to poor fact checking and documentation this kind of demonstration tricks seem to exist. What you basicly do is this: You have an original clear image. You add some artificially generated noise to it and you apply a noise recovery algorithm to remove it! Aaaand ta daaaa perfect recovery is possible, in real life or just for an artifically generated noise?


How do I measure frequency of an audio file into an array? So I know how to plot fft the audio file within Matlab and a window pops up. My question is, what next? I want to know how to go about measuring the frequency(s). I have two similar audio files, of two people saying the same thing, and i want to see if there's a relationship (correlation) between their voices in terms of frequencies and speed (bpm)

Speech processing is a very long and detailed topic, but let's make a start.

Correlation

This probably won't work, but the first thing you can try is cross-correlation of the signals.

xc = xcorr(speaker_1, speaker_2);

Note that really all you may get out of this is the delay between each sample (if you're lucky).

You might also get some level of correlation coefficient by taking the peak value of

xc = xcorr(speaker_1/std(speaker_1), speaker_2/std(speaker2));

Where Next?

If that doesn't work (I'm not hopeful), you will need to look at lots of other options:

Blocking the data into frames
Finding features within these frames (LPCs, etc.).
Using some method to time-align these between speech samples (DTW)
Finding some measure of "likeness" between the time-aligned features.


What is Typical RTP Packet size for H.264 video transmission I am encoding my video using JM reference software. The encoded file is divided into RTP packets. I can set a RTP packet size in bytes using configuration file. But i am unable to find what should be the RTP packet size for H.264 video used in video telephony? I know RTP packet size is variable but there should be some limit.
Any help would be highly appreciated.

RTP Payload Format for H.264 Video

   0                   1                   2                   3
   0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1
  +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
  |V=2|P|X|  CC   |M|     PT      |       sequence number         |
  +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
  |                           timestamp                           |
  +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
  |           synchronization source (SSRC) identifier            |
  +=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+
  |            contributing source (CSRC) identifiers             |
  |                             ....                              |
  +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

    RTP header according to RFC 3550

The RTP header information to be set according to this RTP payload
   format is set as follows:

Marker bit (M): 1 bit
      Set for the very last packet of the access unit indicated by the
      RTP timestamp, in line with the normal use of the M bit in video
      formats, to allow an efficient playout buffer handling.  For
      aggregation packets (STAP and MTAP), the marker bit in the RTP
      header MUST be set to the value that the marker bit of the last
      NAL unit of the aggregation packet would have been if it were
      transported in its own RTP packet.  Decoders MAY use this bit as
      an early indication of the last packet of an access unit, but MUST
      NOT rely on this property.

     Note: Only one M bit is associated with an
     aggregation packet carrying multiple NAL units.  Thus, if a
     gateway has re-packetized an aggregation packet into several
     packets, it cannot reliably set the M bit of those packets.

Payload type (PT): 7 bits
      The assignment of an RTP payload type for this new packet format
      is outside the scope of this document and will not be specified
      here.  The assignment of a payload type has to be performed either
      through the profile used or in a dynamic way.

Sequence number (SN): 16 bits
      Set and used in accordance with RFC 3550.  For the single NALU and
      non-interleaved packetization mode, the sequence number is used to
      determine decoding order for the NALU.

Timestamp: 32 bits
      The RTP timestamp is set to the sampling timestamp of the content.
      A 90 kHz clock rate MUST be used.

H.264 also implements B-Frames. This type of frame fills in information between P-Frames. This means that the B-Frame will need to be held until the next P-frame arrives before the information can be used. B-Frames are not used in all modes of H.264. The encoder decides what type of frame is best suited. There are typically more P-frames than I-frames. Lab analysis has shown TelePresence I-frames to generally be 64K wide (50 packets @ 1316 bytes), while P- frames average 8K wide (9 packets at 900 bytes).

Hope this helps


Is there some literature to understand finite state ISI channel? I am trying to understand channels with memory, where the state of the channel gives the past input information. I want to delve further into capacity for channels with ISI. I came across a presentation given here, on the $5^{th}$ page of it I found a Markov and dicode representation like this:
.

I want to know if there is some primarily literature that concerns itself with understanding this channel, and specifically how to read this binary channel with ISI, and also with dicode channels. Any help will be great

you can also start by "signal processing for everyone" for Gilbert Strang, it's brief and easy to understand. Then you can go with Proakis, "Digital Communications". also there are very good online videos on youtube for India teachers, you are follow them.


Higher side lobe levels in linear array beam pattern I am working on beam forming of high frequency (>400 KHz) linear array for imaging sonar. I have gathered sensor data of 80 channel linear array sub merged halfway in 8 m deep water tank with a transmitting probe 6 m away from it at same depth. I recorded multiple pulses and perform offline analysis on it. The problem i am facing is very high side lobe levels in a beam pattern developed using recorded data. Can any one kindly guide me what could be possible causes of higher side lobe levels?

Thanks & Regards

Nauman

There are several practical items that may cause reduced beamforming performance and lead the higher sidelobes. You haven't said how high exactly these sidelobes are, nor what levels you expected.
Some common areas to investigate further:

Array element spacing. To prevent the presence of grating lobes your element spacing needs to be less than $\frac{\lambda}{2}$. You haven't indicated the bandwidth of your signal, so I'll use the
400 KHz you gave and assume a speed of sound to be 1500 m/s. Then your spacing needs to be less than $\frac{1500}{2\cdot 400\cdot10^3}=1.875$ mm. Of course you should use the highest frequency of signal of interest for this calculation.
I'm assuming that you are using IQ sampling. Phase errors in the IQ demodulation signal or slew in the clocking (synchronization) across the 80 channels may degrade the performance. At this high frequency you will need very small tolerances to limit the effect of these phase errors.
Have you considered the effect of Multipath propagation? The effect of the signals bouncing off the side of the tank or the water surface can lead to signal levels with a higher intensity (constructive interference) than the level of your direct path.


Minimal number of additions in convolution The Winnograd algorithm can be used to reduce the number of multiplications in convolution. Is there a known method of reducing the number of additions in convolution?

One method which attracted my attention recently is Discrete Time Random Sampling.

Which is an approximation method useful to reduce both multiplications and additions in filtering.

Filter and $N$-length signal by a $L$-length filter can be performed in $\mathcal{O}(LM)$  instead of  $\mathcal{O}(LN)$. Where $M\ll N$.
Which is done by picking samples at random in a way described in the paper.

See also Frequency-Shaped Randomized Sampling for a more sophisticated approach which allows you to shape the spectrum of the error.


How to Improve Signal Processing Skills? I've heard that signal processing is fundamental field that can be applied to many other fields (MRI, EEG, audio, speech recognition, communication even to economics)

So, I want to contribute myself to studying this field deeply.

I already took Signals and systems and DSP, and also in this semester, I'm taking Statistical Signal Processing.

This may be stupid question
But, I cannot form any idea of how I can improve and study signal processing.
Taking a theoretical theory is not enough to get the concept.
And even if I learned some basics, I still don't know what signal processing is in general.
Constructing a filter or classifier is main goal of signal processing?

Could you share the general insight and recommend the method I can follow?

One of the best ways is to hang out on this community.
Read questions and answers.
Try to replicate results on answers.

Once you have more knowledge, try answering questions of others.
If you encounter a question on a field you don't know, find resources to read about it until you can answer the question (Assuming it is not deep and requires deep knowledge).

Other classic approaches would be books.
Yet, I recommend Hand On approach.
Namely if you read a book, choose those which has code inside them to replicate results and have better understanding of real world models.

You may even look for some courses on Coursera:

Digital Signal Processing.
Audio Signal Processing for Music Applications.
Fundamentals of Digital Image and Video Processing.

There are few options on edX as well:

Discrete Time Signal Processing.
Signals and Systems, Part 1.
Signals and Systems, Part 2.
Discrete Time Signals and Systems, Part 1: Time Domain.
Discrete Time Signals and Systems, Part 2: Frequency Domain.

For more classic recommendation you may refer to What Resources Are Recommended for an Introduction to Signal Processing (DSP)?

Enjoy the ride!


Why does my I/Q data from radar are not centered? I am working on the IQ raw data from FMCW radar. I plot the I and Q data, they seem to be strongly shifted from the center. And the amplitude of I and Q also varies. I don't understand why it looks like this. I think, they should local on the center and with fixed amplitude.
Thank you very much if someone can help.

The IQ data from an FMCW radar would have the frequency difference from the delay between the transmitted and reflected signal given the transmitted signal is a changing frequency with time; so the result is what should be expected. Given the ramp rate of frequency the round trip delay from target can be determined from the frequency of this result.

The result here appears to have a stronger higher frequency AND lower frequency tone as the result of two different targets. The lower frequency that is superimposed would be due to the closer target.


Why is this matrix invertible in the Kalman gain? In the wikipedia article about Kalman filters, the well-known expression of the matrix of Kalman gains is given:
$$ \mathbf {K} _{k}=\mathbf {P} _{k\mid k-1}\mathbf {H} _{k}^{\text{T}}\mathbf {S} _{k}^{-1} $$
with
$$\mathbf{S}_k=\mathbf {H} _{k}\mathbf {P} _{k\mid k-1}\mathbf {H} _{k}^{\text{T}}+\mathbf {R} _{k}.$$

I understand that $\mathbf{R}_k$, as a covariance matrix, can be asked to be non-singular: it is reasonable to believe that no variance is zero. But this does not answer my question: why is $\mathbf{S}_k$ invertible?

Note that $\mathbf{P} _{k\mid k-1}$, just like $\mathbf{R}_k$, is also a covariance matrix, and for this reason it is (at least) positve semi-definite, i.e., $\mathbf{y}^T\mathbf{P}_{k\mid k-1}\mathbf{y}\ge 0$ for $\mathbf{y}\neq\mathbf{0}$. Now set $\mathbf{y}=\mathbf{H}_k^T\mathbf{x}$ to see that also $\mathbf {H} _{k}\mathbf {P} _{k\mid k-1}\mathbf {H} _{k}^{\text{T}}$ is at least positive semi-definite (positive definite if $\mathbf{P}_{k\mid k-1}$ is positive definite and $\mathbf{H}_k$ has full rank). Finally, note that the sum of two positive semi-definite matrices is positive semi-definite.

For invertibility, we require that the sum of the two matrices is positive definite. This is the case if at least one of the two matrices is positive definite. In practice we can rather safely assume that both $\mathbf{R_k}$ and $\mathbf{P} _{k\mid k-1}$ are positive definite. However, $\mathbf {H} _{k}\mathbf {P} _{k\mid k-1}\mathbf {H} _{k}^{\text{T}}$ is usually only positive semi-definite due to rank-deficiency of $\mathbf{H}_k$, so for invertibility of the sum of the two matrices we have to rely on the positive definiteness of $\mathbf{R}_k$.


Should you include division by coefficient of variation in WV reconstruction for its normed 1D signal? I make a new ECG time series from WV spectrum of original signal and its L2 energy normalisation.
I am thinking if the reconstruction step benefits from covariance at each time point.
I take later a norm of the matrix, which creates the new time series.

result there is no benefit of adding cov postop into the L2 new signal - verified, strong evidence
now, testing if the application of cov on the initial energy spectrum helps

Nonlinearity

I think the coefficient of variation should not be applied to the initial spectrum because there exists a bijection between the two data sets, although a nonlinear one.

  The problem here is that neither sensor is better in this case, because the data sets are direct conversions of each other, but the CV of each data set is different: 0.79 versus 0.42.

Hovever, I am not sure if the nonlineary of the problem case makes the thing different.

Source: https://en.wikipedia.org/wiki/Coefficient_of_variation

There is no benefit of dividing the new improved signal by the covariance of the original ECG and its improved signal like the following code; by the two decimal accuracy in Corrgrams and other more precise tools (eigen cor)

% Covariance http://dsp.stackexchange.com/q/35291/5633
% https://en.wikipedia.org/wiki/Coefficient_of_variation
utp = utp / sum( cov(u,utp) )';

Fig. 1 Left-hand-side figure without division by covariance and Right-hand-side figure with division by covariance - about marginal dependence and independence of the cases

Other inspected methods where I see no difference

Eigen cor matirx of the two
-M^{-1} about conditional independence and dependence

Please, let me know how to get more decimal digits in R corrgram panel.cor; I am trying to replace the problem by Didzis' p-value corrgram discussed there but having problems with finding the upper bound for the amount of points to maintain significant p-value.


is there a mathematical relation between percent overlap and window length or amplitude I am trying to plot short time energy of an utterance with percent overlap as an input variable.

As I referred to the short time energy equation it can be represented as square summation of $x(n)*w(n-m)$ where $w(n)$ seems to be window length.
But I couldn't figure out how to write it as percent overlap.

Percent overlap merely sets how often you measure the energy, relative to the window width.  You can re-measure every sample, or once per annum, or any density in between.


Shannon formula, why square the channel gain? In Shannon formula below:
$$R = W\log_2\left(1 +\frac{P_t\cdot H^2}{W N_0}\right)$$
why square the channel gain? When calculating the channel capacity, why do we need to square the channel gain?

Notation: $W$ is the bandwidth, $P_t$ is the transmit power, $H$ is the channel gain and $N_0$ is the thermal noise power density.

The Shannon channel capacity equation, as I remember is

$$\begin{align}
R &= W \log_2 \left( 1 + \frac{S}{N} \right) \\
\\
  &= W \log_2 \left( \frac{S+N}{N} \right) \\
\\
  &= W \big( \log_2(S+N) \, - \, \log_2(N) \big) \\
\\
  &= \frac{W}{6.02} \big( (S+N)_\text{dB} \, - \, (N)_\text{dB} \big) \\
\end{align}$$

where $S$ is the signal power and $N$ is the noise power (both over the same bandwidth $W$).

Noise power for a given bandwidth (thin enough that noise power density $N_0$ is considered constant is

$$ N = N_0 W $$

and the received signal power is the transmitted power, $P_t$, times the power gain of the channel (assuming both are constant over the bandwidth $W$).  The power gain of the channel is the square of the gain in amplitude because instantaneous power is the square of instantaneous amplitude.  So the received signal power is

$$ S = |H|^2 P_t $$

assuming all of the transmitted power is in the channel bandwidth $W$ and the channel gain, $H$, is constant over the bandwidth $W$.


Hidden Markov Models - Distinct Observation Symbols I really, really hope someone can help me with regards to this question.

I'm trying to implement a Hidden Markov Model (based off this paper: Here

I understand the processes, but, I do not understand what M would represent in the data I am trying to train the HMM with.

I am given this example:

"N = the number of hidden states
 M = the number of distinct observation symbols
 T = the number of observations

 So, for the English text example, if you let N = 2, M = 27 (26 letters
 plus word-space), and T = 50,000 (number of input letters to use), you
 should see that the 2 hidden states correspond to consonants and
 vowels."

This example works for the English Dictionary, I understand that. BUT I am attempting to train the HMM with MFCC Coefficients of a file (stop.mfc) which contains 4k+ values. Now my interpretation would be that: T = 4000; (The size of the 'Observable' sequence) and N = 2; ("Stop" and "Go") so therefore what would M represent in the example I am giving? Which is: Differentiating between someone saying "Stop" or "Go" Would M infer the number of training samples I have?

I really hope someone can help me with this.

Firstly, this is not because you have two words to identify that you need $N = 2$ states. Your goal is not to train a model with two states - one for each word to recognize - but to train 2 models, one for each word to recognize, and each of these models will have as many states as necessary. In fact, each state in your HMM should correspond to a distinct "stage" in the pronunciation of a word - and will very likely correspond to a phoneme. Your vocabulary size (here, two: "stop" and "go") is external to this. For "stop", there are 4 phonemes. For "go", there are 2 phonemes. So you train a 4-state left-to-right model on the "stop" data; and independently of this, a 2-state left-to-right model on the "go" data. To recognize a word given its MFCCs, you evaluate which of these two models has the highest likelihood given the data. If you had to recognize words within a lexicon of 10 words, you would similarly train 10 HMMs, one for each word, each of these models having a number of states suitable to the length/complexity of the word to recognize.

You need to step back and ask yourself "why HMMs in the first place?". We need HMM for speech recognition because words are made of a sequence of distinct elements in sequence (phonemes). If we want to describe/recognize the word "stop", we need to learn a description which is expressive enough to capture that "first it sounds like ssss, for a short while, then it is tttt for a short while, then it is oooo for a longer amount of time, then it is pppp for a short moment". HMMs are a good match for expressing that - states are phonemes, the transition matrix (which will be here diagonal + upper diagonal) indicates that we move through the word from first phoneme to last phoneme, staying a variable amount of time in each phoneme, and the distribution associated with each state indicates how each phoneme translates into your acoustic features.

It seems also that you are mixing up discrete HMM (in which the observations are drawn from a discrete distribution associated with each state) with continuous HMMs (in which the observations are scalars or vectors, characterized by a continuous distribution such as a gaussian). So the parameter $M$, number of distinct observation symbols, is irrelevant in your case, since your observations are 13-dimensional vectors, an uncountable set! ($M$ would be... the cardinality of the continuum).

I am afraid the introduction material you have picked is not directly relevant to speech recognition - though it is useful for applications in which HMMs are used to recover hidden structure from a discrete observations (and there are many of them, for example parsing/tagging in NLP). Try to master this material without thinking much about your speech recognition problem, and then move on to material about continuous HMMs with multivariate normal distributions - and finally to continuous HMM with mixtures of multivariate normal distributions (since this is what is likely to work best for speech).


What is the purpose of the following code? So I'm attempting to setup a LED visualizer using a python library called Impulse. I've already repaired some of the initial issues with the code in my own custom fork, and now I'm attempting to change the chunk size to something a bit more high-resolution. I got the idea from this question, but the issue is I'm a bit lost as to what a certain part of the code does.

I've kinda hit a wall. The bit of code is here. I don't understand what it refers to by fft_max, but I need to figure out what it is so I can determine how to calculate and expand the array to match the higher sample chunk. I apologize if this is in the wrong section, as it is highly related to programming, but I figured someone here might know what fft_max might be referring to.

Here is a plot of fft_max:

I think fft_max is a kind of equalization table so that the visualization would be "flat" over all frequencies for a typical input. It does not look carefully crafted as it is so spiky. Probably the author just used some example audio input and kept book of the largest encountered values of the magnitudes of the bins, hence the name of the table. If you need a different length fft_max, it would be fine to interpolate the original fft_max say by linear interpolation, or to roll your own (not a cigarette).

You could use the synthetic pink power spectral density $a^2/f$, where $a$ is a normalization constant and $f$ is frequency. In terms of magnitude which you need to put into code, that would be $a/\sqrt{f}$, the square root of power. Put some reasonable value to $f = 0$, for example copy the value at the first non-zero frequency.

Figure 1. Pink magnitude spectrum (red diamonds) compared to the original fft_max (blue dots).


Why is the Z-transform so important in digital filters analysis and design? Please elaborate on why this mathematical transform can help analyzing as well as designing any type of digital filter.

The Z Transform is to discrete-time (digital) signals precisely the same role that the Laplace Transform is to continuous-time (analog) signals.

Linear Time-Invariant (LTI) Systems (a.k.a. "filters"), are made up of signal-processing elements that fall into 3 fundamental classes:

adders (devices that add two signals).
scalers (devices that scale a signal by a constant).
"reactive" elements (devices that are able to discriminate w.r.t. frequency).

Element classes 1. and 2. are essentially the same for analog or digital filters.  The are sometimes called "memoryless" devices or elements.

For an analog filter (or "analog LTI system"), those reactive elements would be capacitors or inductors.  They integrate ($s^{-1}$) one signal to become another.  That turns a sine signal into a cosine signal or shifts the phase by $\pm$ 90°.

$$\cos(\Omega t) = \sin(\Omega t + \tfrac{\pi}{2})$$

For digital filters, the reactive elements are delay elements.  A unit delay (a delay of exactly one sample period $T$) will delay any signal, including a sinusoid, by 1 sample or $T$ units of time ($z^{-1}$).  That shifts the phase by an amount that is dependent on frequency

$$\sin(\Omega (t-T) ) = \sin(\Omega t - \Omega T)$$

or

$$\sin(\omega (n-1) ) = \sin(\omega n - \omega )$$

Any LTI system that acts as a "filter", a device to filter out some frequency components and leave others, must have reactive elements (or "non-memoryless" elements or components having memory) in order to discriminate one frequency from another.  And such a filter will shift phase which will normally be different for different frequencies.  But a memoryless LTI system (which is just a scaler) will not discriminate between frequencies nor will shift phase, except for possibly by 180°, which is just a polarity reversal or scaling by a negative constant.


Step-by-Step calculation of stationary wavelet transform gives shifted results? Good morning,

I am doing an SWT by hand to understand it better, and have a couple of questions. Maybe someone here has experience with this?

I am using the book "Conceptual Wavelets" by Fugal, where they also perform an SWT by hand. This is the page I'm referring to:

Now the details coefficients we get in the book are
-80 0 0 0 80 0 0 0 0.

When I did the convolution by hand I got this same result. Matlab agrees also.

However, when I do a level 1 SWT in Python or Matlab, the details coefficients look like this:
0 0 0 56.6 0 0 0 -56.6

I understand that the numbers are different - the book uses [1 -1] as wavelet filter, and most algorithms use [1/(sqrt(2)) -1/(sqrt(2))]. (I don't know why, it's not explained in the book.)

But the Matlab result is shifted to the left! Somehow, Matlab gets rid of the delay in the result during the computation of the details coefficients. They also reduce the convolution back to the right length. In the book, all of these transformations are done at the end, to get the original data back.

So my questions are:

Is the delay always 1, for every 1D SWT?
How does Matlab calculate its result, and how can I replicate this in my calculation? Is it always a shift to the left? Which element do I drop to reduce the size of the result? A simple convolution seems to be not enough.
(And if someone happens to know a source that mentions why haar filter values are so often 1/(sqrt(2)), I'd be happy to learn about that, too.)

Thanks in advance

Why normalized filters

Last question first, "Why $1/\sqrt2$": Because it makes the (Euclidean) norm of the filter $1$, so that the whole wavelet filter bank operation (if done right, that is, the decimated version) is orthogonal/isometric. It is a design choice, there is nothing wrong in staying with the integer values and correct the combined factor during reconstruction.

The book example is not really the standard wavelet transform

However, the book you cite is doing the wavelet transform undecimated. Which is a valid mathematical transform, but leaving out the decimation step, the down- and up-sampling, will lead to data bloat. By applying the two filters, the amount of data representing the signal is doubled, so some of that is redundant. This redundancy can be reduced by transmitting (or keeping for further processing) only every second element of these sequences.

Undecimated or otherwise redundant transforms give freedom in how to represent the signal, that is, in some non-linear post-processing one can set some of the data to zero and compensate (more-or-less) lossless-ly by changing other coefficients. However here I think this example is just a teaching aid.

Decimated version of the examples

So taking of cD1 and cA1 only every second entry gives cD1(2:2:end)=[0 0 0 0] and cA1(2:2:end)=[160 160 0 0], together they encode the same signal, which now has more zeros and is thus better suited for compression. In the reconstruction, zeros will be added to fill the gaps, and then the convolution with the reversed/biorthogonal filters is applied. So

[1 -1]*[0 0 0 0 0 0 0 0] = [0 0 0 0 0 0 0 0 0]

and

[1 1]*[160 0 160 0 0 0 0 0] = [160 160 160 160 0 0 0 0 0]

so that in the sum the result is twice the original signal, $(\sqrt2)^2$, with a zero added at the end.

If one shifts the whole process by one place, the reduced signals are [ -80 0 80 0] and [80 160 80 0] so that in the reconstruction one has to add

[1 -1]*[-80 0 0 0 80 0 0 0] = [-80 80 0 0 80 -80 0 0 0]

and

[1 1]*[80 0 160 0 80 0 0 0] = [80 80 160 160 80 80 0 0 0]

to give [0 160 160 160 160 0 0 0 0]. Here we get, additional to the factor 2, also the shift, or an additional zero at the start of the sequence. However, by using the first, third etc. elements in down-sampling, we used implicitly an extension of the sequence S to the left by zero. This additional zero element is what is obtained now explicitly in the reconstruction.

In general wavelet transforms, the wavelet filters have lengths larger than 2, leading to more "interesting" boundary effects that include a shift by about the length of the filter.


Gabor filter bank for gist features I was reading a paper about Gist features. I was interested about its code implementation which i found it here). Into LMgist.m there is a function called createGabor, this is the code:

function G = createGabor(or, n)
%
% G = createGabor(numberOfOrientationsPerScale, n);
%
% Precomputes filter transfer functions. All computations are done on the
% Fourier domain.
%
% If you call this function without output arguments it will show the
% tiling of the Fourier domain.
%
% Input
%     numberOfOrientationsPerScale = vector that contains the number of
%                                orientations at each scale (from HF to BF)
%     n = imagesize = [nrows ncols]
%
% output
%     G = transfer functions for a jet of gabor filters

Nscales = length(or);
Nfilters = sum(or);

if length(n) == 1
    n = [n(1) n(1)];
end

l=0;
for i=1:Nscales
    for j=1:or(i)
        l=l+1;
        param(l,:)=[.35 .3/(1.85^(i-1)) 16*or(i)^2/32^2 pi/(or(i))*(j-1)];
    end
end

% Frequencies:
%[fx, fy] = meshgrid(-n/2:n/2-1);
[fx, fy] = meshgrid(-n(2)/2:n(2)/2-1, -n(1)/2:n(1)/2-1);
fr = fftshift(sqrt(fx.^2+fy.^2));
t = fftshift(angle(fx+sqrt(-1)*fy));

% Transfer functions:
G=zeros([n(1) n(2) Nfilters]);
for i=1:Nfilters
    tr=t+param(i,4);
    tr=tr+2*pi*(tr<-pi)-2*pi*(tr>pi);

    G(:,:,i)=exp(-10*param(i,1)*(fr/n(2)/param(i,2)-1).^2-2*param(i,3)*pi*tr.^2);
end

if nargout == 0
    figure
    for i=1:Nfilters
        contour(fx, fy, fftshift(G(:,:,i)),[1 .7 .6],'r');
        hold on
    end
    axis('on')
    axis('equal')
    axis([-n(2)/2 n(2)/2 -n(1)/2 n(1)/2])
    axis('ij')
    xlabel('f_x (cycles per image)')
    ylabel('f_y (cycles per image)')
    grid on
end

If you run this function with the input parameters or = [8 8 8 8] and n = 192 (as in the original implementation) it should create this image:

This image shows the contours of the gabor filters in the frequency domain for z = 0.6, z= 0.7.

Unfortunately this gabor filter implementation is different from what i've seen so far. My questions are:

Look at the matrix param, the second column definitely seems to be
the frequency scale, where fmax = 0.3 and c = 1.85 is the constant.
The fourth column definitely seems to be the orientation scale. The first column is a constant of value 0.35; the only reference i've found in that paper about 0.35 is this sentence (page 8, right column):

  The model of Eq. (7) provides correct fitting for all the eight
  categories for frequencies below 0.35 cycles/pixel (as noise and
  aliasing corrupt higher spatial frequencies, see Fig. 4)

So what the first and third column represent?
Which brings me to this line of code:
G(:,:,i)=exp(-10*param(i,1)*(fr/n(2)/param(i,2)-1).^2-2*param(i,3)* pi *tr.^2);
this implementation is very different from what i've seen so far, i
don't think i recognize the formulation. In particular the second
term of the exponential: 2*param(i,3)* pi *tr.^2 controls the
orientation scale and also the particular "stretch" of the gaussian
in the frequency domain as shown in the picture, which I think is
used to get a better filter overlap in the frequency domain. Can
somebody clear up this implementation?

Here is my understanding of these points. I apologize if my answer comes too late. I post it anyway in the hope it helps someone.

The 4 parameters in the param vector respectively controls the radial width, radial location, angular width, and angular location of the filter in the Fourier domain. In other words, the filter's response is a (separable) Gaussian function of the frequency modulus and orientation (that is, the polar coordinates). The parameters #2 and #4 gives the location of this Gaussian while parameters #1 and #3 gives its scale. However, I have no explanation regarding the "magic" values (maybe they were simply obtained by trial & error).
These filters are not Gabor filters since their Fourier transforms are not Gaussian functions of the 2D frequencies. On the figure you showed, Gabor filters would result in ellipsoid curves instead of these drop-shaped curves. I guess the authors made this choice to have a better coverage of the spectral plane.


Filter response to sinusoidal input signal This question is related to Does 46 dB gain of a filter for a frequency imply 200 times more amplitude? In  the accepted answer it is mentioned that "For a sine wave input the amplitude of the steady state output is simply the amplitude of the input multiplied with the magnitude of the transfer function at that frequency"

Above relation is derived as an integral; as a convolution between filter and the input. But in practice we have all discrete signal, therefore can we get 200 times more in practice?

The fact that a sinusoidal signal's amplitude is multiplied by the magnitude of the system's frequency response (evaluated at the input frequency) is a general property of LTI systems, regardless of whether they are continuous-time or discrete-time systems.

This is a consequence of the complex exponential being an eigenfunction of LTI systems, and the value of the frequency response evaluated at the input frequency is the corresponding eigenvalue.

For more details on eigenfunctions of LTI system check out this question and its answers. And I'm sure that there's even more information here on this site.


Evaluating the autocorrelation function of a telegraph signal to get exponential decay? I've been driving myself absolutely crazy trying to do autocorrelation analysis on some data that I've got.  I've read from multiple sources that the ACF of a random telegraph signal / random telegraph noise is given by an exponentially decaying function where the decay constant is the average rate of switching between states (cf. Example 7.3.2 here for example).  Since I had trouble getting anything that even resembled exponential decay with real data, I decided to generate a fake signal that telegraphs between the value 0 and 1 with lifetimes centered binomially around two characteristic values for the 0 and 1 state.  The autocorrelation function for this signal is still not at all exponential.  In each case, the ACF ends up strongly oscillating.

I'm working in Mathematica, using CorrelationFunction to get the autocorrelation function.  Based on my reading, this is a standard definition for autocorrelation.

So my question is can someone show me an example with a numerically obtained ACF from a random telegraph signal that decays exponentially?  Am I missing something stupid here?  Am I using an incorrect definition for the autocorrelation function?

A semi-random telegraph signal is one for which $X(0) = 0$ and $X(t)$ equals $0$ or $1$
according as there have been an even number or an odd number of arrivals in a Poisson Process in the interval $(0,T]$.  A simple physical model would be the least significant
bit (LSB) of a digital counter that is counting the Poisson arrivals.  A random telegraph signal is almost the same except that $X(0)$ is equally likely to have value $0$ or $1$.
(Feel free to translate this to signals having values $\pm 1$ if that is more convenient
or more in accordance with your notation).

With this as background, you are encountering problems because your simulation is
not modeling the random telegraph signal correctly. You should instead simulate a
Poisson process and obtain the random telegraph signal from that.


IQ generation separation To generate I & Q samples from data samples S
I = S(1:2:length(S)-1) * cos(2*pi*f*t) ;
Q = S(2:2:length(S)  ) * sin(2*pi*f*t) ;

I & Q to integer arithmetic
I = I * (2^16) ;
Q = Q * (2^16) ;

I Q combined
y = bitshiftleft(Q, 16) ;
y = bitor(y, I);

and mixed with carrier fc
y = y * cos(2*pi*fc*t) ;

To recover the I & Q, we do:
after mixing y with the carrier fc to get baseband y, we do
Q = bitand(y, 0xFFFF0000) ;
Q = bitshiftright(Q, 16) ;
I = bitand(y, 0x0000FFFF) ;

I am assuming dividing is the reverse of the above multiplication
I = I / cos(2*pi*f*t) ;
Q = Q / sin(2*pi*f*t) ;

My question is the above  correct especially dividing by cos and sin?
I am not sure how to calculate t in all the above.
Please help.

Your approach is putting the I and Q samples on a "sub-carrier" at frequency $f$. If that is really what you wanted to do, the typical demodulation approach is to multiply by the sub-carrier again and low pass filter to get I or Q. Dividing by cosine or sine creates time samples that go to infinity since you will divide by zero with those functions. The multiplication operation produces a baseband component and a component at twice the carrier frequency; this is made clearer with a quick refresh of the following trigonometric identities:
\begin{align}
\cos(\alpha)\cos(\beta) &= \frac{1}{2}\cos(\alpha-\beta)+\frac{1}{2}\cos(\alpha+\beta)\\
\sin(\alpha)\sin(\beta) &= \frac{1}{2}\cos(\alpha-\beta)-\frac{1}{2}\cos(\alpha+\beta)\\
\sin(\alpha)\cos(\beta) &= \frac{1}{2}\sin(\alpha-\beta)+\frac{1}{2}\sin(\alpha+\beta)
\end{align}
So in your case:
$$
\cos(2\pi ft)\cos(2\pi ft) = \frac{1}{2}\cos(0)+\frac{1}{2}\cos(2\pi ft) = \frac{1}{2} +\frac{1}{2}\cos(4\pi ft)
$$
The $\frac{1}{2}$ represents the amplitude of I and the $\frac{1}{2}\cos(4\pi ft)$ is the doubled frequency component that is removed with a low pass filter.
However I am not sure that you intended or need to have such a "sub-carrier"? If you were going to simulate QPSK modulation for example, you can modulate directly to your final carrier as in:
\begin{align}
I_t &= I \cos(2\pi f_c t), \quad\text{where}\quad I = {1,-1}\\
Q_t &= Q \sin(2\pi f_c t), \quad\text{where}\quad Q = {1,-1}
\end{align}
And further the I and Q samples are sent simultaneously by adding the modulated outputs $I_t$ and $Q_t$ above (hence the increased data rate; we send two bits in each "symbol).
To demodulate, assuming we have done the carrier recovery work to have a good estimate of $f_c$ with all phase offsets removed, we can multiply the received sequence by $2\cos(2\pi f_c t)$ and low pass filter to recover I, and multiply the received sequence by $2\sin(2\pi f_c t)$ and low pass filter to recover Q:
Noting the demodulation steps below where LPF() denotes removing the $2f_c$ frequency component with a low pass filter:
\begin{align}
I &= LPF\big\{2 I_t\cos(2\pi f_c t)\big\}\\
Q &= LPF\big\{2 Q_t\sin(2\pi f_c t)\big\}
\end{align}
The Q component goes to zero on the I output since the result is in the $\sin(\alpha)\cos(\beta)$ form resulting in $\sin(0) = 0$:
$$
0 = LPF\big\{2 Q_t\cos(2\pi f_c t)\big\}
$$
And similarly the I component goes to zero on the Q output:
$$
0 = LPF\big\{2 I_t\sin(2\pi f_c t)\big\}
$$
Therefore we can "multiplex" I and Q together given this ability to separate the two that operate simultaneously and independently in time.
Note I highly prefer to think and model using exponential frequency terms instead of sines and cosines - it is so much simpler! In fact this best explains implementations with I and Q paths as well, since to implement a complex signal you will need two real sequences $(e^{j\omega t} = \cos(\omega t) + j\sin(\omega t))$. Please see Frequency shifting of a quadrature mixed signal for more details on that.


Applying periodogram-based PSD For plotting a Power Spectral Density graph based on FFT results,
I had used two different approaches, but they produced very different results, and it'd be great if anyone knows why.

Note: From this point on, $N$ is assumed to be the total number of data samples, $f_s$ is assumed to be the sampling frequency.

APPROACH #1:

(from here: www.mathworks.com/help/signal/ref/periodogram.html    )
$$PSD(f) = \frac{\Delta t}{N} \left| \sum_{n = 0}^{n = N-1} x_n e^{-i2 \pi fn}\right|^2$$

of which I assume the $\Delta t$ to be the same as $1/f_s$, with it being multipled by 2 for most cases (one sided periodogram). So the idea was implemented as:
$$PSD(f) \cong \frac{2}{f_s N} \left| FFT(x_n)  \right|^2 $$

Note: As per one of the responses, APPROACH #1 would've been more accurately interpreted as
$$PSD(f) \cong \frac{1}{f_s} \left| FFT(x_n)  \right|^2 $$

APPROACH #2:

(from here: http://www.bitweenie.com/listings/power-spectral-density-matlab/)
$$PSD(f) = 20 \log_{10} \left( \frac{|FFT(x_n)|}{\sqrt{N f_s}}\right) + 30$$

This approach was implemented with no interpretative modifications.

RESULTS:

Using approach #1, my results looked like a typical FFT, while approach #2 produced the expected log-scale result.

QUESTION:

However, the approach #1 from MATLAB seems like the more commonly-accepted way to do the periodogram PSD.
Is there a fault somewhere in the interpretation of this approach to cause the discrepancy in results?

It's exactly the same method, only the second is converted to (power) dB ($\text{dB}:\,x\mapsto 10\log_{10}x$):

Let's convert your first formula to dB:

$$\begin{align}
P_{1,\text{dB}} &= 10\log_{10}P_1\\
&=10\log_{10}\left(\frac{\Delta t}{N}\lvert\text{FFT}\rvert^2\right)\\
&=10\log_{10}\left(\frac{\lvert\text{FFT}\rvert^2}{f_sN}\right)\\
&=10\log_{10}\left(\frac{\lvert\text{FFT}\rvert^2}{\sqrt{f_sN}^2}\right)\\
&=10\log_{10}\left(\left(\frac{\lvert\text{FFT}\rvert}{\sqrt{f_sN}}\right)^2\right)\\
&=10\cdot2\cdot\log_{10}\left(\frac{\lvert\text{FFT}\rvert}{\sqrt{f_sN}}\right)\\
&=20\log_{10}\left(\frac{\lvert\text{FFT}\rvert}{\sqrt{f_sN}}\right)\\
\end{align}$$

The $+30$ just happen because your seconds source is considering $\text{dBm}$ (decibel of milliwatt) instead of $\text{dBW}$ (decibel of watt), assuming your input values are in Volt over a impedance of 1Ω or so.

Note: you nowhere in your question stated the physical units of what your input data is; be very careful when "implying" something can be represented as $\text{dBm}$ or $\text{dBW}$. That makes a statement about your signals origin (namely, that your ADC or wherever the signal came from is calibrated, or will be converted to analog with calibrated equipment) if you're dealing with real-world observations/signals. This assumption is usually wrong – if you just got a good microphone and a good sound card (or a good receiver and a good SDR device) your digital amplitudes will be proportional, but not identical, to the voltages observed, or even more significantly, to the sound intensity (or EM field strength) that caused the electrical signal in the first place.


motion compensation for range migration algorithm Is there a simple to follow explanation of how to do second order motion compensation for range migration algorithm?

The one described in Carrara et al is for the ideal case. For the Extended Omega K algorithm, how does one relate the deviation in platform position to the amount of phase compensation to apply in terms of kx?
Is it a scaling by $ 2\pi/\delta x $ where $ \delta x $ is the deviation along the azimuth, and the slant range deviation is assumed to have been compensated away via a mix of phase compensation and/or interpolation and/or pixel shifts?

I have been trying to incorporate MoCo to the extended Omega-k as well following Reigbar's paper (https://ieeexplore.ieee.org/document/1648550). I understand that for the first order MoCo, the range deviation needs to be calculated from a reference range (i.e. image center) and for the second order, the differences between the deviation for all ranges and the deviation for the reference range need to be used. However, calculating the deviation in the above recommended method didn't work for me. Especially, I'm not sure if I calculated the deviation for the second order MoCo right. I used r = (c*fast_time)/2 for calculating all the range positions and then subtracted those from the non-ideal platform positions in order to know their range deviations. Is that correct?

Thank you.


Recursive systematic convolutional code (RSC) realization in MATLAB How can I implement $\rm RSC(2,1,2)$ in MATLAB? Can I use convenc function to build RSC?

The convenc function has 2 arguments: code = convenc(msg,trellis) where trellis=poly2trellis(ConstraintLength,CodeGenerator), and the generator polynomial for this RSC is:
$$
G=\left(1,\frac{1+D^{2}}{1+D+D^{3}}\right)\quad\text{or}\quad G=\left(1,\frac{5}{7}\right)
$$
but for poly2trellis argument CodeGenerator need to be in octal.

The question is already answered but let me structure it for future readers.
From the MATLAB help for the function poly2trellis, here is how you can do it:

TRELLIS = poly2trellis(CONSTRAINTLENGTH, CODEGENERATOR, FEEDBACKCONNECTION)

Hence for the figure in the question:
 RSC_trellis = poly2trellis( 3, [7 5], 7);

Or, for the other possible RSC encoder, i.e. choosing the lower branch to be fedback:
 RSC_trellis = poly2trellis( 3, [7 5], 5);


Discretizing a continuous time signal This may be a very straightforward question but I've been trying to understand this for a long time and can't wrap my head around it.

How would one go about discretizing the following continuous time signal:

signal_continuous = 1/6 * sin(880*pi*t) + 1/15*sin(1760*pi*t);

As far as I understand, if I wanted to sample the above signal with a sampling rate of 8Khz over two seconds. I would do something like this (in Matlab speak):

t = [0:1/8000:2-1/8000];
signal_discrete = signal_continuous./t;

Is this the whole picture or have I completely missed the point here?

My rationale behind doing the above is that since the continuous time signal is, well, continuous and if I want to sample it at 8000 cycles per second over a two second period. I would just divide the signal by that (t in this case). However, what happens to the rest of the signal? Or rather, does this not distort the whole signal (if the original was longer than 2 seconds, for example)?

EDIT:

I was complicating things for no reason. In Matlab I would just need to do the following:

t = [0:1/8000:2-1/8000];
discrete_signal = 1/6 * sin(880*pi*t) + 1/15*sin(1760*pi*t);

This would create a two second discrete-time signal out of the original continuous signal.

The theory provided in the answer, I think, is beneficial so I'm leaving the question up.

A continuous-time signal $x(t)$ is sampled to create an ideal discrete time signal $x[n]$ by just taking the instantaneous value of the continuous-time signal at the sampling instants. If you sample the signal at an interval $T_s$, then:

$$
x[n] = x(nT_s) = [0, x(T_s), x(2T_s),\ \ldots \ ]
$$

In your case, $T_s = \frac{1}{8000}\text{sec}$. In practice, there are no ideal samplers, so this model doesn't strictly hold, but it is close enough for most applications and is the model that you would typically see used in introductory signal processing textbooks and courses.