What is the difference between continuous, discrete, analog and digital signal? It's my first time studying DSP and I've faced a problem finding a convenient definition.

Are the following definitions correct? And if so why there are some resources defining it in other terms such as "Digital signal: is a signal with discrete time and discrete amplitude"

Discrete time signal: X-axis (time) is discrete and Y-axis (amplitude) may be continuous or discrete.
Continuous time signal: X-axis (time) is continuous and Y-axis (amplitude) may be continuous or discrete.
Digital signal: Y-axis (amplitude) is discrete and X-axis (time) may be continuous or discrete.
Analog signal: Y-axis (amplitude) is continuous and X-axis (time) may be continuous or discrete.

A signal is indeed a function. Given a signal $f(x)$, according to whether continuous or discrete for both the variable $x$ and the function $f(x)$, there are four types of combinations:

(1) $\mathbf{continuous}$ $x$ and $\mathbf{continuous}$ $f(x)$

This is the most common $\mathbf{analog}$ signal.

(2) $\mathbf{continuous}$ $x$ and $\mathbf{discrete}$ $f(x)$

For this one, we can imagine the ideal base-band waveform used in digital communication. As this one:

(3) $\mathbf{discrete}$ $x$ and $\mathbf{continuous}$ $f(x)$

This is indeed the signal in most the "digital signal processing" textbooks. An example, as others have pointed out, is the output of the CCD sensor.

(4) $\mathbf{discrete}$ $x$ and $\mathbf{discrete}$ $f(x)$

This is the $\mathbf{digital}$ signal. Digital signals are used in practical implementation aspects, and they actually exist in a conceptual manner.

If we concern the discrete feature of the function, the problem will be more complex, therefore, in most "digital signal processing" textbooks, the signals are $\mathbf{not}$ digital indeed. An interesting fact is that, for the the classical textbook by A. V. Oppenheim, the name was "digital signal processing" in the 1st edition, but the name was changed to "discrete-time signal processing" for the later editions.


Why delete high frequence in fft return the area of the most changes? I'm new to ||image-processing|| and || fft filters ||

I have a relative smooth $\mathtt{2D}$ (circle) gaussian surface (height $10$) with several "holes" in it, and the surrounding area [has] some noise (height below $0.5$).

Here's what I did, I transfer[ed] the matrix $\mathcal{A}$ with $\mathcal{A}_f= \text{fftshift}( \text{fft}(\mathcal{A}))$, then I used the for loop to delete all the high frequency component, delete the component with magnitude compar[able] to DC component ($\approx 2.5\mathtt{e}5$), i.e. if $\mathcal{A}_f(i,j) > 1\mathtt{e}3$ , then [set] $A_f(i,j)=0$.

Then, I did a backwards transformation $\mathcal{A}_\text{result}=\text{ifft}(\text{ifftshift}(\mathcal{A}_f))$.

What's so interesting is that, $\mathcal{A}_\text{result}$ clearly identified all the holes (the region with the most changes) as [spikes], and the rest of the graph is almost $0$ with some noises.

My question is that:

What happened here? Why the image returns the value of holes as [spikes]?
Is there any other way to identify the holes with Fourier or convolution filter?

This line of yours

  delete the component with magnatude compare to DC component(=2.5e5),
  i.e. if Af(i,j)>1e3 , then Af(i,j)=0

indicates an effective highpass filter on the typical smooth image data.

Indeed for typical smooth images high, frequency coefficients are generally (but not always) smaller in magnitude compared to the low frequeny coefficients. So deleting large magnitude coefficients will effectively delete low frequency coefficients. But this is not exactly same thing with deleting low frequency coefficients directly. neverthless, such an operation would return mostly a darkened image with certain edges and contours outlined, essentially a very a primitive edge detection exemplified with the following code:

I = im2double(imread('Lena.bmp'));
figure,imshow(I);

Ik = fft2(I);
th = 100;
Ik( abs(Ik) > th) = 0;
figure,imshow( 3*real(ifft2(Ik)));

with the result:


Which Approach Is Better for Decomposing an Image into High Frequency and Low Frequency Components? Which approach is better or there is mathematical justification for using Bilater filter and Fourier Transform to decompose a image into High Frequency and Low Frequency Component.

Both Bilateral Filter and Fourier Transform can be used for getting High and Low Frequency component of an image.

Is one approach is better than the other or are both same in getting high and low frequency component of a image, or are they used in some particular context of image.

I know that Bilateral Filter uses range weights and Spatial weights into consideration for getting high and low frequency components.

The bilateral filter is just that: a filter. It does not decompose the image into anything. It just attenuates "medium" frequency content in a clever manner. In other words, it smooths things that are mostly smooth further, and leaves sharp edges sharp. What you are left with has more low and high frequency content (on a relative basis) than what you started with.

The Fourier transform, on the other hand, actually does reveal the entire spectrum for an image. Unlike a filter, it does not change the information content of an image. It just rearranges it into the frequency, rather than spatial, domain. The bilateral filter does not reveal the spectral content of an image, so really only the Fourier transform (or one of its related transforms, like cosine and wavelet) can be used to get the high and low frequency components of an image.


Formulas to plot Frequency Impulse Response Coefficients I made this question https://dsp.stackexchange.com/questions/19731/newbieplotting-fft-and-impulse-response-coefficients-using-java

I have the IRC: Impulse Response Coefficients:

$h(k), k = 0, 1,.., N-1$

Low Pass: $h_D(k)=\sin(k\omega_c)/{k\pi} , k\neq0; h_D(0)=2f_c$

High Pass: $h_D(k)=-\sin(k\omega_c)/{k\pi}, k\neq0; h_D(0)=1-2f_c$

Band Pass: $h_D(k)=(\sin(k\omega_{cH})-\sin(k\omega_{cL}))/{k\pi}, k\neq0; h_D(0)=2(f_{cH}-f_{cL}) $

Band Stop: $h_D(k)=(\sin(k\omega_{cL})-\sin(k\omega_{cH}))/{k\pi}, k\neq0; h_D(0)=1-2(f_{cH}-f_{cL}) $

Hamming Window Function:

$w(k)=0.54+0.46\cos( \frac{2\pi k}{N} )$

so we have:
$h(k) = w(k)\times h_D(k)$

obviusly my $h(k)$ can be represented as:

double[] h = new double[N];  //N:odd integer
...
//hD[(N-1)/2] -> center or hD(0)
...
for (int  k = 0; k < N; k ++) {
  h[k] = hD[k]*w[k];
}

I want to get the formula to evaluate the h[k]

In matlab to do this:

[hFlt,fFlt] = freqz(h,1,256);
mLP = 20*log10(abs(hFlt));
plot(fFlt,mLP);

I want to obtain hLP and fLP values to plot filters response (without using matlab).

But I don't know the formulas(procedure or algorythm) to plot the frequency response.

Assuming you have a set of (windowed) filter coefficients $h[n]$, $n=0,1,\ldots,N-1$, and you want to evaluate the corresponding frequency response at $ M\ge N$ equidistant normalized frequency points

$$f_k=\frac{k}{2M},\quad k=0,1,\ldots,M-1\tag{1}$$

then you need to evaluate the following formula, which is the discrete Fourier transform (DFT):

$$H[k]=\sum_{n=0}^{N-1}h[n]e^{-j\pi nk /M},\quad k=0,1,\ldots, M-1\tag{2}$$

The actual frequencies in Hz are given by $f_k\cdot f_s$, where $f_s$ is the sampling frequency.

Equation (2) is exactly what you get with the Matlab command

[H,w] = freqz(h,1,M)

where $f_k=$ w(k+1)/(2*pi) (because Matlab indices start at $1$, not at $0$).

If you want to evaluate the frequency response at arbitrary, not necessarily equidistant normalized frequencies $f_k$, you need to evaluate

$$H[f_k]=\sum_{n=0}^{N-1}h[n]e^{-j2\pi nf_k},\quad 0\le f_k\le 0.5\tag{3}$$

which for the choice of $f_k$ according to (1) of course reduces to Eq. (2).


Fourier transform of moving image sensor Consider an image sensor that is moving with some velocity $\mathbf{v}(t) \in \mathbb{R^3}$, and an orientation $\mathbf{R}(t) \in SO(3)$. Treat those vectors as fixed. At a specific time $t$ the projection of the scene on the sensor is:
$$u_{\mathbf{v},\mathbf{R}}(t)$$
The sensor pixel data is now sampled in a series fashion. This means that each pixel will be sampled at a different time $t$.

Now suppose that you want to calculate the Fourier transform of a frame. How would you do it? Can you do it without interpolation? I'm after accuracy here.

Furthermore, what is the effect of the lack of sample synchronization on the FFT i.e. If you just take the FFT without first synchronizing the pixels. Is there a rigorous mathematical treatment of this subject? What is the frequency representation of the sampled data?

Edit

This is for a rectangular image sensor in an interferometry application. This means that although the velocity is small $v \approx 1 \mathrm{m/s}$, the image changes rapidly.

I thought about it some more and I guess what you would have is aliasing for each pixel equal to the frame rate, and each pixel will have a different amount of phase shift equal to the pixel clock.

I feel that the question needs a little bit more clarifying but there probably is enough material here to provide some sort of response and adjust it later if needed.

To cut a long story short: Rolling Shutter distortion. This concept is key here (irrespectively of whether or not you have a rolling or global shutter camera that is unsynchronised to the interferometer).

Now, I will treat $v, R$ as independent from time here for a minute. Also, $F_{FR}$ will stand for frame rate.

Your "problem" seems to be that due to $v$, each pixel row might actually be sampling spectral lines from two radically different wavelengths. If this is true then each one of your frames will acquire a slant whose angle will be proportional to $\frac{v}{F_{FR}}$.

In the time dependent version, the frames are capturing a "curve" whose curvature is proportional to $\frac{v(t)}{F_{FR}}$ and...good luck with that.

Since you are dealing with interferometry, you are only really imaging ONE single focal plane, therefore, you could apply a simple (anti)-skew transform to counteract the skew caused by the rolling shutter. This, of course, is within limits. In any case, the rolling shutter effect of the camera inserts a phase error in the image which, to an extent, can be counteracted but in the end it all comes down to controlling $v$.

Assuming that your camera can maintain a constant $F_{FR}$ for the duration of the scan, then essentially, you don't have a 2D camera.

I would like us to imagine it as one long linear $M \times N$ camera (because of $v$ and interferometry) only it is wrapped around with a "wrap" factor that depends on $v$.

So, theoretically, you could find a $v$ which scans a spectrum range in one frame. Maybe $v$ is such that the spectral range is contained in two frames or three frames or four...Obviously, the more, the better. Therefore, each pixel row now has a $\lambda_{start}$ and a $\lambda_{end}$. So, I want you to imagine a chirp that is wrapped across the camera sensor or in fact across many different frames.

Therefore, you could convert your single $M \times N$ image, into one long vector and then do a spectrogram on this vector. This vector is essentially encoding wavelength to space (across the camera sensor).

Does this have any distortions? Of course it has, but it would be like trying to obtain the frequency response of a system by feeding it a chirp. Within limits it is perfectly possible. If you scan it too fast, you might skip a spectral line

Perhaps it is easier to see why $v(t)$ would insert a phase error to the "chirp" when considering the camera as one long image sensor (?).

Hope this helps.


Fixed point multiplication with negative numbers I am stumped on a simple problem. Let's say I have two 4 bit numbers in Q0.3 format. One sign bit and three fractional bits. So I can represent $-1$ through to $0.875$.

Let's now say I wish to do this calculation: $-0.25 \times 0.875$. Which is:

$$
\frac{-2}{2^3} \times \frac{7}{2^3}
$$

Which means I am multiplying $1110$ ($-2$) by $0111$ ($7$). Of course the answer is $-0.21875$ or $-0.25$ using the closest Q0.3 number.

Let's do the working.

$$
1110 \times 0111 = 01100010
$$

which when viewed as a Q0.6 number is $1.100010$, which is $-0.46875$ by my books. Why is this incorrect? I expect an answer of $1.110010$ ($-0.21875$).

What have I done wrong?

When multiplying two's complement numbers, you have to perform sign extensions to the operands to meet the number of digits your multiplication will yield, i.e., in your case $4 + 4 = 8$ digits.

$$
11111110_2 \times 00000111_2 = 11110010_2
$$

As there are $2 * 3$ fractional bits, the result is $1.110010_2 = \frac{-14}{2^6} = -0.21875$. Normalizing this number to $3$ fractional bits in Q0.3 format yields $1.110_2 = -0.25$.


Mel filter in MFCC - is it necessary? Is it necessary to use filter bank in MFCC process? Can anyone explain what is the Mel filter? I know that frequency in hertz is converted into Mel scale but is this formula can be directly applied after the Fourier transformation of the speech signal? Why the filter bank is needed?

Mel filter bank is important due to following reasons:

It applies the Mel-frequency scaling, which is perceptual scale that helps to simulate the way human ear works. It corresponds to better resolution at low frequencies and less at high.
Using the triangular filter-bank helps to capture the energy at each critical band and gives a rough approximation of the spectrum shape, as well as smooths the harmonic structure. In theory you could manipulate on raw DFT bins, but then you are not reducing the dimensionality of your features - this is the whole point of doing filter-bank analysis, to capture the spectral envelope.

Obviously warping of frequency scale can be altered. Other mappings are possible, such as Bark, linear, etc. Mel is just one out of many options. To finalize, you can always calculate more filter-banks to have coarser scale (getting closer to actual DFT), but you will notice that at some point it will degrade the performance and your features do not carry essential information.


Loop bandwidth for symbol timing recovery I am trying to design the Symbol timing recovery (STR) block of DVBS Receiver.

How should I chose the loop bandwidth and damping ratio for determining the loop filter coefficients $k_1$ and $k_2$?

In some literature they are given as $B_nT_s= 0.005$, $\zeta = \frac{1}{\sqrt 2}$, while in others it is given as $B_nT_s= 0.05$, $\zeta = 1$.

How should I proceed for choosing the $B_n$ and $\zeta$?
Should $B_nT_s$ factor change as we change $T_s$ (say from $0.1$ to $10^{-5}$ sec)?

There is an optimum loop BW that maximizes the SNR, and this applies to both Symbol Timing Recovery as you inquire about as well as Carrier Recovery. The specific answer depends on the characteristics of the noise source involved in your system, such as clock jitter and phase noise, and the modulation you are using; specifically how the signal energy is distributed in frequency. Without knowing those an exact answer cannot be determined, but I can at least share some insights from my own experience with it.

To state generally, I have found that the optimum loop BW to usually be somewhere between R/100 and R/20, where R is the symbol clock rate being recovered.  I start with a value in this range and then through simulation with a valid noise and waveform model (or lab measurement) determine the actual optimum using SNR metrics vs Loop BW setting.

As for damping factor, I also start with a damping factor of 0.5 and then through evaluation using a step or impulse response test adjust for possible ringing due to parasitic delays that may be present in the loop.

In addition to the SNR considerations given above, there are also acquisition time constraints or possible dynamic noise conditions that would also motivate you to have as fast of a timing loop as possible: The faster the loop, the faster the acquisition time, and recovery time from an interference event, or the ability to track a highly dynamic change if it exists.  There may also be waveform specific constraints such as the duration of a header in a packet based waveform used to establish timing. However as the BW for the tracking loops increase beyond a significant fraction of the timing symbol clock rate (such as >R/20), tracking errors will start to dominate the noise in SNR, since the loop will start to track out the modulation of interest.

To explain further from the perspective of a trade of Loop BW and SNR, the Loop BW sets the boundary at which your recovery loop will track timing noise (jitter/phase noise) in your signal.  Timing noise has components at all frequencies (such as that viewed with a phase noise spectral density), and the lower frequency components of the timing noise (the change in zero crossing locations that is moving very slowly) will be tracked by the timing loop, and therefore suppressed. However the higher frequency components will not be tracked by the timing loop and therefore remain and contribute to noise as part of SNR. The timing loop can be viewed as a high pass filter on the jitter/phase noise components, passing the higher frequency components (including the modulation signals) and supressing the lower frequencies at a rate depending on the order of the loop. Thus you see the motivation to make the timing loop as fast as possible in the interest of tracking out as much timing jitter as possible. However, some of the signal energy of interest will be part of the "noise" being supressed by the timing loop; and this can be quantified by understanding the spectral characteristics of the modulation used. As the loop BW increases, too much signal energy is also removed, and thus it is clear how there can be an optimum setting for the Loop BW that maximizes SNR.

This hopefully provides insight into the approach and what to look for in the development of the system. In conclusion, I set the Loop BW as fast as possible (typically in the range of R/100 to R/20 where R is the symbol rate) while monitoring through simulation with accurate noise and waveform models the SNR metric; and either choose the point of best SNR, or if other factors such as acquisition time or more stringent dynamic conditions are involved, maximize BW without going below a target SNR metric.


Reconstruction of a signal using an ideal DAC I think I was able to do part (a) correctly, but now am stuck on part (b). I don't know how to approach this part or what is the theory behind it.
My MATLAB code is shown below. Please let me know if I have any mistakes.

% part a
T = 0.2;
t = 0:0.001:3*T;
xc = 10 + 3*sin(20*pi*t+pi/3)+ 5*cos(40*pi*t);
% for sampling time t = 0.01n
n = 0:0.01:3*T;
xn = 10 + 3*sin(20*pi*n+pi/3)+ 5*cos(40*pi*n);
figure(1)
plot(t,xc);
hold on
stem(n,xn,'filled');

% part b
% constructing yr to produce x[n]

This is an ill posed question.

First it doesn't talk about the output DAC reconstruction sampling rate which is necessary for the reconstruction filter design and which also determines the frequency spectrum scaling of the reconstructed signal. Anyway then you may assume the DAC sampling-rate to be equal to ADC sampling rate for simplicity, which are both $F_s = 100$ Hz in your case. (btw your matlab code is wrong here, you have taken the sampling period to be $0.001$ which should be $0.01$.)

Furthermore, it talks about an ideal DAC which cannot be realized by any code, but, its output can be mathematically described perfectly without any computation by definition of the ideal DAC to be $$y_r(t)=x_c(t)$$ as long as there is no aliasing during the ADC stage which is the case here. So it's not clear if they want you to actually compute an approximate output $y_r(t)$ or simply evaluate the exact described output $x_c(t)$ in terms of input samples, which is a null computation.


Question regarding estimation of signal tone parameters (frequency, amplitude, and phase) using Macleods algorithm I plan to use Macleods algorithm in order to estimate the parameters (frequency, amplitude, and phase) of one or more sinusoids from a block of N uniformly spaced samples.

Before applying the mentioned algorithm, I would like to know how to pre- or postprocess the captured signal. Can you recommend the following two actions in order to ensure that Macleods algorithm delivers best results?

Capture more data points of the stationary signal to increase the resolution of the estimated parameters?
Postprocess the captured signal with a (hanning) window before
applying the algorithm

Thank you for your answers.

I'm not sure about Macleod, but generally windowing is not good for that style of estimator. See this paper for details about similar estimators.

Generally the more data, the better for all parameters. Provided you use an estimator without bias, it should do better with more data.

It depends on what  you mean by "resolution".

If you mean the standard deviation of the estimate, then that will be governed (at higher SNRs) by the Cramer-Rao lower bound for frequency estimate $\hat{f}_0$, which is:
$$
{\rm var}(\hat{f}_0) \ge \frac{12}{(2\pi)^2 \eta N(N^2-1)}
$$
where the signal model is
$$
x[n] = A \cos(2\pi f_0 n + \phi) + w[n]
$$
the amplitude is $A$ and phase $\phi$ (which are also both unknown), the number of samples is $N$, $\eta = A^2/(2\sigma^2)$ is the SNR and $w[n]$ is white, Gaussian noise of variance $\sigma^2$.

See equations (3.41) of Kay's Fundamentals of Statistical Signal Processing, Volume 1: Estimation Theory.


How to implement Fourier Descriptor of an image? I want to implement Fourier Descriptor of an object. I have read link. However, I have some questions about normalizing Fourier Descriptor.

First, if I want to normalize the position of the starting point, according to the above link, This is done by subtracting the phase of the second Fourier descriptor $\phi_1$ from the phase of all Fourier descriptors and weighting by $k$; that is:

$a(k)=a(k)e^{-j \phi_1 k}$

Is the weight $k$ equal to the $k$ in $a(k)$ ?

Second, how could I do in order to achieve rotational invariance when I implement Fourier descriptor ? Thank you very much for your reply.

Yes, these k are the same.

In the basic sense, the rotational invariance is obtained by only taking the magnitude of the descriptors. However, in practice this might lead to information loss and decreases discrimination power. For this reason it is typically prefered to augment the phase information with additional terms. There are many ways to do this but for a simple approach follow:

Invariant Fourier Descriptors Representation of Medieval Byzantine Neume Notation

Dimo Dimov, Lasko Laskov


Determine if two signals are the same? I'm pretty new to signal processing. My data is the $(x,y,z)$ acceleration of someone walking with a smartphone in his/her pocket.

Goal

Given two signals, determine if it's from the same person twice, or if it's from two separate people.
Have some idea of the degree of confidence in the decision from 1.

Attempt

Segmented the data into the cyclic portions (each step and retraction of foot), and acyclic portions (maybe waiting between steps or when the user stops to talk to someone).
Calculate the statistics of the cyclic portions, and acyclic portions (mean, s.dev, median, cross correlation, time between each step and such).

Having said that, I lack the base knowledge for determining if two signals are 'similar' to one another which I'd imagine comes from some distance function. I've found ResearchGate thread but I was wondering if there's anything else that I should include (especially since my goals are not the same as that users). The methods they mention are:

cross-correlation
covariance
mutual-information
autocorrelation

My background is mostly stats / ML / software so I'd really appreciate it if any explanations were dumbed down!

Valuable results can be obtained using Mel Frequency Cepstral Coefficients (MFCC) and Dynamic Time Warping (DTW).

Divide the signal in multiple segments. I recommend to use an overlap of 50%. Compute the MFCCs for each frame, then you will obtain a vector of MFCCs for each frame. Finally you will have a matrix of [No of MFFCs, Number of frames], for each signal. To compare the similarity between sequences with different time length, compute the Euclidean distance between the two matrices. Also you can use the standardized Euclidean distance (check scipy.spatial.distance.cdist for Python).

cost_matrix = cdist(x, y, metric='euclidean')

Thus, you will obtain a cost matrix. Applying the DTW algorithm on this cost matrix, you will find the score which will tell you the similarity between the input signals. I recommend you to see how a different number of MFCCs affects the results and also take in consideration the length of the signals.

I recommend this link to read more details about MFCCs and there is a useful video tutorial about DTW. Also you can find many articles about DTW, it is a well known measure for the similarity between two time series.


What is the transfer function of this block diagram I and my friend have different answer for this block diagram

Mine:
$$y[n] = -\frac23x[n] + x[n-1] -\frac12[n-2]$$

Hence

$$H(z) = \frac{1}{-\frac23 + z^{-1} -\frac12z^{-2}}$$

My friend:
$$q[n]=x[n]-\frac12q[n-1]$$
$$y[n] = -\frac23x[n]+q[n-1]$$
$$Y(z) = \frac23X(z)+z^{-1}Q(z)=\frac23X(z)+z^{-1}\frac{X(z)}{1+\frac12z^{-1}}$$

Hence

$$H(z)=\frac{Y(z)}{X(z)}=\frac23+\frac{z^{-1}}{1+\frac12z^{-1}} = \frac{4z^{-1}+2}{3+\frac32z^{-1}}$$

I don't know how can he come to that result but I'm not sure about mine either

As pointed out in the other answers, your friend's approach is correct. However, his end result for $H(z)$ isn't.

With $q[n]$ the input signal of the delay element, you get $q[n-1]$ at its output, and the given difference equation for $q[n]$ follows easily. The output $y[n]$ can then be written in terms of $x[n]$ and $q[n-1]$. In the $\mathcal{Z}$-transform domain you get algebraic equations in terms of $X(z)$, $Q(z)$, and $Y(z)$, from which you can eliminate $Q(z)$, leaving you with an expression for $H(z)=Y(z)/X(z)$:

$$H(z)=\frac{Y(z)}{X(z)}=-\frac23+\frac{z^{-1}}{1+\frac12 z^{-1}}=-\frac23\frac{1-z^{-1}}{1+\frac12 z^{-1}}\tag{1}$$

So it was basically just a sign error that messed up the end result.


What is Magnitude and Phase actually represents in Fast Fourier Transform? When we reconstruct the image from its magnitude only using Inverse Fast Fourier Transform, then why the resulting image looks nothing like but the original image? Can someone please explain me in detail, what is the role of magnitude and phase in FFT and IFFT?

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

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

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


Estimation of Two Closely Spaced Frequencies? What is the best frequency estimation algorithm for two closely spaced frequencies in term of the minimum frequency spacing achieved?

This answer is a consideration of the frequency resolution when using the DFT but does not answer the question specifically in terms of what is the best algorithm. (see comments by Amro).

If the two closely spaced signals are similar in magnitude (and sufficiently stronger than any other signals, i.e. high SNR), then the highest frequency resolution (which is what you are interested in) is achieved with a rectangular window. This means simply take an FFT with no windowing, and the signal bandwidth will be 1 bin wide or similarly equal to 1/T where T is the time length of your composite signal. This relationship holds regardless of what algorithm is used; the frequency resolution is related to the time length of the signal with the best resolution of 1/T. If the two closely spaced signals are not similar in magnitude, you will run into dynamic range issues as the rectangular window has the worst sidelobe levels, which can mask lower level signals. In this case we use windowing; specifically we multiply the time series by a window function prior to taking the FFT, which will significantly reduce sidelobe levels, but will have the drawback of increasing the width of the main lobe (and hence decrease frequency resolution).

See this answer for further related info: Specific frequency resolution


Should I observe a signal multiple times for MUSIC algorithm? In a paper Multiple emitter location and signal parameter estimation, the well-known MUSIC algorithm starts by calculating covariance (autocorrelation, whatever) matrix with:

where overlines stand for expectation operation, i.e., $\overline{XX^H}\triangleq E\left[XX^H\right]$. As I understand, the MUSIC algorithm needs to observe signals multiple times to get $S$, i.e., $S=\frac{1}{N}\sum_{i}^{N}{X_{i}X_{i}^{H}}$
On the other hand, when I googled implementations of the MUSIC algorithm, most of them does not perform expectation operation. I mean, they just observe signals only one time and perform the remaining process of the algorithm, i.e., $S=XX^H$
Is there a reasonable explanation that the MUSIC algorithm does not need to observe signals multiple times, but only need to observe signals only once?

This question is very much related to this answer.

The issue is that the expectation operator to get the auto-correlation (auto-covariance) is generally replaced in signal processing with the sample auto-correlation.


Adjustable Notch Filter help I'm converting the simulink diagram found in this paper to C code.  I'm not familiar with matlab/simulink(too expensive) so need help interpreting the diagram.  The problem I experience is that once a value like $\theta$, or $x_1$ becomes zero, zero is propagated everywhere due to the multiplication.  The paper states that the integrators are initialized to zero, so I don't see how the system works.

Any ideas?

Below is the simulink diagram and my C code:

The block diagram is found in the paper above as Figure 1 and Figure 2 :

typedef struct  { float zeta,gamma;} PARAM_ANFfloat;
typedef struct  { float err,regressor,x1,theta;} STATE_ANFfloat;

typedef struct {PARAM_ANFfloat p;STATE_ANFfloat s;} ANFfloat;

void dspANFfloat_perform(ANFfloat *anf, float y)
{
    anf->s.err = y - anf->s.regressor;

    anf->s.theta = -anf->p.gamma * anf->s.x1 * anf->s.theta * anf->s.err;

    anf->s.x1 = anf->s.theta * anf->s.err * 2 * anf->p.zeta - anf->s.x1 * - anf->s.theta * anf->s.theta;

    anf->s.regressor = anf->s.x1;

}

In this line:

  anf->s.theta = -anf->p.gamma * anf->s.x1 * anf->s.theta * anf->s.err;

you seem to be conflating $\theta$ and $\dot{\theta}$ (the derivative).

Equation 5 of the paper is:

$$
\dot{\theta} = -\gamma x_1 \theta \left[ y(t) - \sum_{l=1}^n \dot{x}_l \right]
$$

Note that the left-hand term is the derivative of $\theta$, not $\theta$ itself.


Audio noise detection with python Have many audio files 1 minute long. Some of them are normal. Some of them are noise.

Here is a normal file:

This is a noise:

And this is noise too:

How do auto detect noise using python?

BTW: if you suggest any external console program which can do that it will be cool too.

First, there are many ways to do your task. Depending on your background, some might be easier than others for you to understand and use.

If you are familiar with signal processing, you might construct a noise model and based on this noise model to detect whether or not a given signal is of pure noise. This type of method often involves hypothesis tests, and test statistics. Based on these things, it rejects or fails to reject a given sample accordingly. This is a quite traditional way to solve your problem.

If you really donot care about details here, you might complete ignore signal and audio. Purely consider what you have are different data streams, and you can use all kind of machine learning methods to differentiate two classes of things, in your case, one is noise and the other is non-noise. If you are going to adopt this approach, python sklearn package will be your good friend with many great tutorials teaching 0-level people to do machine learning.

Finally, let me remind you that one big difference of the later approach from the former one is that you need prepare labels of each signal sample. Otherwise, ML algorithms will have a better chance to differentiate these two classes when you have labels and many samples. Anyhow, both approaches are belonging to the pattern recognition family, and thus they are definitely related to each other. For example, you might view a probabilistic noise model is somewhat a naive bayes classifier.


Finding n Amplitudes by DFT, what is correct normalization Please forgive me if this has already been asked.
Let us assume an example with $x(t) = \sum_{i=1}^N A_i \sin(2\pi f_i t) $ under a given sampling frequency $f_s$, frequencies $\omega_i$ and Amplitudes $A_i$. My ideas was to calculate the DFT transform $\tilde{x}(\omega)$, find the peaks of $|\tilde{x}(\omega)|$ and infer the original amplitudes from them.

However, there are several obstacles.

There is leakage, and the Maximum of the peak is a bad indicator for the Amplitude of the signal; the normalization factor seems to be dependent on several factors.
The Peaks overlap.
I need a window function to get rid of the side lobes, which again changes my normalization.

Do you have experience in doing this? What is the best practice? Any useful literature?

Do you want to get the peak, or the amplitude, of the time domain signal given the spectrum? Since energy is conserved, the time domain energy amplitude can easily be found by adding the energy from each frequency bin of the spectrum. However, you cannot find the time domain peak simply by looking at the FFT frequency bin magnitudes. Consider an extreme case: the spectrum of an impulse is fairly smooth in spite of the fact that the impulse will have a strong peak.


What do the Contour Levels Produced by Matlab Function Contourc Represent I have been trying to understand what the contour levels are.

As far as I understand a contour of the function $Z(x,y)$ is a set of points in the $(x,y)$ plane, such that $Z(x,y)$ is fixed at some constant value, that constant value is the contour level.
Matlab docs state that the contourc function takes an argument that determines the levels upon which the contours are computed.
I wanted to try that out, so I constructed a picture with $10$ rectangles, each of which has a distinct depth starting from $0.1$ ending up at $1.0$, and I ran the contour function to see the highlighted contours according to the specified levels in the input, but the results weren't expected.

For instance, I ran the code on $0.1$ contour level and the contour function produces all the contours in the picture, not only the ones with $0.1$ level.

Here is the Matlab code used to construct the picture and find its contours:

pic = zeros(500, 500);
val = 0.1;
shift = 40;

for j = 1:10
  for i= 50:450
    pic(shift*j,i) = val;
  end
  for i= 50:450
    pic(shift*j+25,i) = val;
  end
  for i = shift*j: shift*j+25
    pic(i, 50) = val;
    pic(i,450) = val;
  end
  val = val + 0.1;

end

imshow(pic);
figure(2);
contour(pic, [0.1 0.1]);

I'd be so grateful if someone could clue me in.

So you:

  I constructed a picture with 10 rectangles, each of which has a distinct depth starting from 0.1 ending up at 1.0

and then:

  I ran the code on 0.1 contour level

The matlab function is doing the right thing. A rectangle of height 0.1 has a contour at 0.1.  A rectangle of height 1 also has a contour at 0.1 (and 0.05 and 0.5, and 0.9999999, but not at 1.000001).  The rectangle exists from 0 to 1 (in this case).


Conversion of airflow signal into 1-D as a input for CNN Trying to detect apnea from airflow signal. I have 100 files of the patient each file contains at least 8 hours of data and sampling rate is 32 Hz. How to prepare my dataset which is used as an input channel. I am new in this domain please help me.

The airflow signal is usually respiratory input in a polysonogram. If your input files are .EDF files, then you need to extract the airflow signal only, which is a 1-D signal. I think there are many implementations of edfread.m that are available online.

https://www.mathworks.com/matlabcentral/fileexchange/31900-edfread


Inter-point interpolation using FIR filter Suppose I have discrete noisy signal $X = (0.096, -0.0632, 0.351, 0.531, 0.360, 0.006, -0.320)$ sampled at discrete time points $T = (1, 2, 3, 4, 5, 6, 7)$. Filtering (zero-padded) $X$ with symmetrical FIR $G = (\frac{1}{4}, \frac{1}{2}, \frac{1}{4})$ gives me a good approximation of underlying $f(t)$ but only at discrete time points $T$. Is it possible to use filter coefficients to approximate the values of $f(t)$ between discrete values, e.g. estimate $f(\sqrt 2)$?

In Interpolation with an FIR filter user Hilmar mentions Whittaker–Shannon interpolation formula, but it is sensitive to noise. Is it possible to construct smooth but more noise-robust interpolation by using FIR coefficients somehow?

Or, separate interpolation step (e.g. spline or polynomial interpolation) after FIR application is the only practical way to go?

From Hilmar's answer to the linked question, the Whittaker–Shannon interpolation formula is:

  $$x(t)=\sum_{n=-\infty }^{\infty }x[n]\cdot \operatorname{sinc}\left(\frac{t-nT}{T}\right).$$

If the signal is low-pass and the noise is high-pass, to filter high-pass noise, we can make a low-pass version of the interpolation formula:

$$x(t)=\sum_{n=-\infty }^{\infty }x[n]\cdot 2f_c/f_s\cdot \operatorname{sinc}\left(\frac{t-2f_c/f_s\cdot nT}{T}\right),$$

where $f_c$ is the cutoff frequency and $f_s$ is the sampling frequency. The highest signal-to-noise ratio is obtained by setting the cutoff frequency to a frequency for which the signal and the noise have equal power. Actually the Whittaker–Shannon interpolation formula is low-pass as well, with $f_c=f_s/2$.


finding the audio pulse compression peak in extra-low SNR scenario? I just want to obtain some general ideas for such a scenario: if the environmental noise is extremely strong and the audio signal at the receiver is extremely weak(for instance, low range propagation), so even after pulse compression(cross correlation using linear chirp signal), the correlation peak is hidden in the noise. Are there any ways to detect the hidden peak?

So this answer is coming from experience with radar systems which use complex signals (in-phase and quadrature signals so that phase information can be recovered), be aware that you may have to adapt it to fit your exact needs:

With radar systems, this is a pretty common problem. Sometimes, for a variety of reasons, the target that is present can’t rise above the noise, even with regular pulse compression gain. In this case, most systems turn to some sort of multi-pulse solution; a common one is called coherent integration. Basically, you send out and receive the exact same pulse N many times.

After doing this, you’ll have N many versions of the same target. Assuming no target migration has happened, you simply add all of the complex signals together. The idea here is that the target signal will cohere across all of them, such that a processing gain is achieved, but the noise will not cohere with itself and it’s power level will remain fairly constant.


How can i know the full duplex does work from the Shannon's theorem? There is a full duplex system below,two user are $A$ and $B$,according the definition of full duplex ,$A$ and $B$ have to transmit and receive the signal at the same time.that is, transmitted signal will become a noise for  the received signal,and there are some formulas below to show the SNR

Let $t_A$ is the transmitted signal from $A$,and $t_B$ is the transmitted signal from $B$.

Now,$A$ receive the $t_B$ from $B$ and transmit $t_A$ to $A$.So the SNR for $A$ to receive the $t_B$ is

$\text{SNR}_{A \text{ decode }t_B}=\frac{\text{power of } t_B }{\text{power of } t_A \ + \text{power of channel noise} + \text{decode noise}}\label{SNR}\tag 1$

And according to the wiki: https://en.wikipedia.org/wiki/Shannon–Hartley_theorem

How do i make sure that $A$ can do decode the $t_B$?

I know the $C_A=\log_2(1+\text{SNR}_{A \text{ decode } t_B})$, but in most of paper, they all also call $C_A$ as a rate, so I think $C_A$ will not be equal to $R_A$.

So how can I know $C_A \ge R_A$,so the $A$ can really receive and decode the  $t_B$?

You think the Shannon-Hartley Theorem is in conflict with Full Duplex (in the sense you describe), and luckily it's not!

The problem really is that your $\eqref{SNR}$ formula isn't relevant for the Shannon-Hartley Theorem (SHT). With respect to the SHT, your $t_A$ isn't noise!

That's easy to explain. SHT makes a statement about the channel capacity $C$, i.e. the maximum mutual information between $B$ and $A$, given the best possible source symbol probability distribution:

$$\label{C} C= \max_{P(B)} I(B;A)$$

Mutual information between the information source $B$ and the information sink $A$ is, in the end, nothing but the statement

  How much uncertainty about $B$ can I remove by observing $A$?

Now, even from reading that, you'll notice that in a perfect transceiver $A$, you simply don't add any uncertainty by adding your own transmission – you know what you're sending, and thus could simply subtract it from what you observe!

Formally, it's easy to show that

\begin{align}
I(B;A) &= H(B) - H(B|A)\\
&=H(t_B) - H(t_B|r_A)\\
&=H(t_B) - H(t_B|t_B+Z+t_A)\text.\label{crucial}\tag{a}
\end{align}

Now, that $\eqref{crucial}$ line is crucial – what's written there literally means

  The amount of information you get about the transmit signal $t_B$ by observing the received signal $r_a=t_B+Z+t_A$ is exactly the expected amount of information (==entropy) of $t_B$, minus the remaining uncertainty about $t_B$ given that you've observed $r_B$, which is the sum of A's own signal $t_A$, independent noise $Z$ and our signal of interest, $B$'s $t_B$.

And that conditional entropy is what the influence of your noise it. Let's look at it a bit closer, $r_A$ be the receive signal at $A$

$$r_A=t_B + Z  + t_A,$$

with $Z$ being your "true" noise.

And hence, since these three are independent:

$$H(r_A) = H(t_B) + H(Z) + H(t_A)\text.$$

Little problem: We're at the end $A$. So, there's nothing uncertain about $t_A$ (we're the one sending that!) and the formula reduces to

$$H(r_A)_A = H(t_B) + H(Z)\text.\tag{ö}\label{known}$$

(Read $\eqref{known}$ as "The entropy of the received signal $r_A$, being $A$, i.e. knowing $t_A$ already").

Automatically that means that $\eqref{crucial}$ also simplifies when we, as $A$,  know our own transmit signal $t_A$:

\begin{align}
I(B;A)_A
&= H(t_B) - H(t_B|t_B+Z+t_A,t_A)\\
&= H(t_B) - H(t_B|t_B + Z)\label{read}\tag{r}\\
&= H(t_B) - H(Z)\\
&= H(B) - H(Z)&\rule{0.5em}{0.5em}
\end{align}

The trick is that in line $\eqref{read}$, the second term is "the amount of uncertainty about one element of a sum, given the value of the sum", and that simply reduces to the amount of uncertainty in the other element.

When you look at that, you'll notice that your Full-Duplex mutual information formula is the same as your (normal) simplex mutual information formula. Hence, the capacity is the same.

In other words: For Shannon-Hartley, any operation that you can mathematically reverse as the receiver, like adding a signal you know, makes absolutely no difference for the channel capacity. Yay!

The fact that there's nevertheless a lot of things to calculate for the Full Duplex use case is because adding a strong signal to a received signal is typically not fully reversible, because the two amplitudes physically don't add up exactly, due to imperfections of the real world (receive amplifier saturation, changed quantization noise, feedback path uncertainties…). But: Your simple SNR mode $\eqref{SNR}$ can't describe that, so it's neither suitable for the ideal case I described above, nor suitable for the non-ideal case. So, if you want to understand in detail why Full Duplex works, and what the limitations are: You'll need to learn quite a bit about information theory!


Sinusoidal Synthesis of audio I am working on a sinusoidal encoder & decoder to compress .wav files. For now I am trying to regenerate the audio by combining all the frequencies of the original .wav file.
Here is the code I am using:

[x fs] = wavread('clip.wav');
x  = mean(x,2);
NFFT = 1024;
hop = NFFT/2;
win = hann(NFFT);

[S F]= spectrogram(x, win, hop, NFFT, fs);
freq = linspace(0, fs/2, hop + 1);
numFrames = size(S,2);
mag = abs(S);
phase = angle(S);

%Decoding
output = zeros(size(x));
t = 0:(1/fs):(((NFFT-1)/fs));
sumOfCosine = zeros(1, NFFT );
frameStart = 1;
frameEnd = NFFT;

%f = frame
for f = 1:numFrames
    sumOfCosine = zeros(1, NFFT );

    %k = sinusoids in each frame
    for k = 1:hop+1
       sumOfCosine =  sumOfCosine + (mag(k,f)*cos(2*pi*F(k)*t + phase(k,f)));
    end

 sumOfCosine = sumOfCosine.*win';
 a = hop;
 % Normalize the window
 sumOfCosine = (sumOfCosine./a);
 sumOfCosine = sumOfCosine';
 output(frameStart:frameEnd) = ((output(frameStart:frameEnd)+ sumOfCosine));
 frameStart = frameStart + (hop);
 frameEnd = frameEnd + (hop);

end

When I play play output, there is still some noise in it. I tried comparing the original matrix x multiplied by the window, win to the output frame by frame, and I can see that output represented by green doesn't completely overlap with the original windowed signal. I am pasting the images for the first three frames below.

As it can be seen there is some error due to the overlapping of windows maybe? I'd appreciate help in fixing this and make the output sound cleaner. Thanks!

For the phases of sinusoids within the frame overlaps to sum constructively, either t has to be adjusted to be frame relative, or the fft phase result has to be adjusted to be relative to the same absolute time as t.


Is it right to talk about dynamic range on frequency domain? I'm comparing three frequency response's plots from simulations and experiments and I observed that two of them look like clamped to a DC component, where the max dB value is i.e. 100 and the min is 40 dB, while on the third one the max is i.e. 120 dB and the min is 10 dB. I want to use the correct term to talk about these intervals.

I've used the term dynamic range to express an amplitude variation interval on time domain, but is it right to use it on frequency domain too?

the dynamic range is the difference between the largest and the smalles value of a quantity can asume, so yes, you can have a dynamic range for the magnitud of your frequency components, I'm not aware if there's another terminology when dealing with the frequency domain. One thing to consider is not only the dynamic range of your signals but also the dynamic range of the variables you'rr using to store it and to calculate the frequency response, altough in the frequency domain (and specially when working with dB) it's not a big issue.


BER calculation in DS-UWB system I need to calculate BER (bit error rate) for a system using spreading spectrum DS over ultra wideband (UWB). Modulation is bpsk. The spreading sequence lengths are not equal for all users. Receiver is matched filter. I know the equations of DS-UWB transmission signal and template signal which I should use at receiver to extract the user's data. I know that I should use convolution between the received signal (transmitted + AWGN noise) and template signal. My problem is, how to calculate the result of that convolution? is it integration or multiplication?
thanks in advance

Assume two discrete-time sequences: $x$ (length $L_\mathrm{x}$) and $y$ (length $L_\mathrm{y}$). Furthermore, $L_\mathrm{x} \geq L_\mathrm{y}$. To calculate the discrete-time (linear) convolution of $x$ and $y$, first pad $y$ with $L_\mathrm{x} - L_\mathrm{y}$ zeros, i.e.
$$
y'(n) =
  \begin{cases}
    y & 0 \leq n < L_\mathrm{y}\\
    0 & L_\mathrm{y} < n \leq L_\mathrm{x} - 1
  \end{cases}
$$
Then calculate the convolution of $y'$ and $x$ as described here, for example.


Validity of an argument that two transfer functions are minimum-phase based on their ratio being minimum-phase Update I think the essence of my question below is this: If the ratio of two transfer functions may be represented exactly as a minimum-phase filter (MPF) plus a pure delay (in the title, I simply refer to this representation as “minimum-phase” for brevity, even though it’s not accurate), does that imply that each transfer function in the ratio must also be represented exactly as an MPF plus a pure delay?
I was reading a paper where an argument is made in the Appendix regarding the validity of representing a transfer function as a minimum-phase filter (MPF) plus a pure delay.
In particular, the goal is to show that two transfer functions $H_L$ and $H_R$ can each be represented as a minimum-phase filter (MPF) plus a pure delay.
That is, if we have $H_L = |H_L| \text{e}^{\text{i}\phi_L}$ and $H_R = |H_R| \text{e}^{\text{i}\phi_R}$, where $\phi_L$ and $\phi_R$ are their respective phase responses, then given that these are rational transfer functions that can be decomposed into a cascade of minimum-phase and all-pass transfer functions, we must show that $\phi_L = \phi_L^{\text{mp}} + \omega\,t_L$ and $\phi_R = \phi_R^{\text{mp}} + \omega\,t_R$, for some constants $t_L$ and $t_R$ (i.e., the all-pass phase response is purely linear), and where the superscript "mp" stands for minimum-phase.
As I understand it, the argument by the authors of the paper goes as follows:
Give $H_L$ and $H_R$ as above, we have $$\ln \left(\frac{H_R}{H_L}\right) = \underbrace{\left(n_R-n_L\right)}_{\Delta n} + \text{i}\underbrace{\left(\phi_R-\phi_L\right)}_{\Delta \phi},$$
where $n_R = \ln |H_R|$ and $n_L = \ln |H_L|$.
Since, in general, we have $\phi_L = \phi_L^{\text{mp}} + \phi_L^{\text{ap}}$ and $\phi_R = \phi_R^{\text{mp}} + \phi_R^{\text{ap}}$, where the superscripts "mp" and "ap" stand for "minimum-phase" and "all-pass", respectively, this means $\Delta \phi = \Delta \phi^{\text{mp}} + \Delta \phi^{\text{ap}} = \mathcal{H}\left(\Delta n\right) + \Delta \phi^{\text{ap}}$, where $\mathcal{H}$ denotes the Hilbert transform operator.
The authors show that this may be expressed as $\Delta \phi = \mathcal{H}\left(\Delta n\right) + \omega t_0$ for a constant, $t_0$.
Let us assume that this equation (which corresponds to Eq. (A8) in their paper) is valid.
The authors then go on to suggest (incorrectly in my opinion) that since $\Delta \phi^{\text{ap}} = \omega t_0$, $H_L$ and $H_R$ must each always be representable exactly as an MPF plus pure delay (when this is, in fact, not known).
I think this argument is flawed because it assumes that $\Delta \phi^{\text{ap}} = \phi_R^{\text{ap}}-\phi_L^{\text{ap}} = \omega t_0$ implies that $\phi_R^{\text{ap}}$ and $\phi_L^{\text{ap}}$ are also each linear functions of $\omega$ (since only then do $H_L$ and $H_R$ each have MPF plus delay representations).
However, this is not the only possibility correct? I think a function like $\omega t_0$, which is linear in $\omega$ can, in general, be obtained as the difference of two non-linear functions of $\omega$, each corresponding to either $\phi_R^{\text{ap}}$ or $\phi_L^{\text{ap}}$, which would mean that an MPF plus delay representation for $H_L$ and $H_R$ as described earlier is not always necessary, right? For example, if $\phi_R^{\text{ap}} = \alpha \omega^2$ and $\phi_L^{\text{ap}} = \alpha \omega^2 + t_0 \omega$, for some scalar $\alpha$, then $\Delta \phi^{\text{ap}} = \omega t_0$. Am I going wrong somewhere?
Thanks in advance for any inputs!

two transfer functions HL and HR can each be represented as a minimum-phase filter (MPF) plus a pure delay.

That is generally not true and it's easy enough to disprove it by counter example. Let's simply look at a first order all pass
$$H(z) = \frac{p-1/z}{1-p/z} = H_{MP}(z) \cdot H_{AP}(z)$$ which has a pole at $p$ and a zero at $1/p$. The "minimum phase" part, $H_{MP}(z)$ is simply unity, i.e. $H_{MP}(z)=1$ and the all pass part is just the original all pass. The all pass phase response is NOT that a of a pure  delay, we clearly have $\Phi_{MP}(\omega) \neq \omega \cdot t_0$
I think you misunderstand the point in the paper: I think they are just saying that a specific class of transfer functions can be approximated well with this model. In this case this class is "HRTF ratios for human heads" and they show that this works particularly well for spherical heads. Since binaural perception is mainly based on interaural differences, approximating the frequency dependent interaural time delay by a single constant value can be quite useful.


What are characteristics of different MR images? MR tools are used in many sectors such as medicine and mining. MR stands for magnetic resonance and MRI for magnetic resonance imaging. What are characteristics of different MR images? What are typical noise to signal ratios and other important decision factors?

Speaking from a medical imaging perspective, the SNR of MR images mainly depends on signal acquisition schemes, image reconstruction parameters, and the MR scanner hardware, such as, number of receive coils. To know about more about other parameters go here: http://wiki.ci.uchicago.edu/pub/HNL/DifficultQuestion/howdoesmriwork.pdf.

Parallel MR imaging is one of the acquisition technique which improves SNR. For example, I'm sharing with you the results from a paper published in journal of magnetic resonance imaging on "Measurement of Signal-to-Noise Ratios in MR
Images" by Olaf Dietrich and others below:

PS: Take the reciprocal to get noise to signal ratio.


How to extract text from a scanned document? The task I am facing is extraction of rotated text (any angle) from a scanned image(attached Image). The image contains a blue sticker on it and need means in extracting the text (printed text) from it. The pattern of the sticker will remain the same for all the time.

As you can see the orientation is not stable and i would be really grateful if anyone could suggest a manner in which i could extract the text.

Reference to the question can be found here, or the technique How to remove text region in a document image? or extract image from a document image

Thanks

I think your workflow should be something like this:

1) Detect the rectangle.

2) Find the rectangle's orientation.

3) Rotate the rectangle and place it in a new image.

4) Feed the rectangle twice to an OCR
(one time with the calculated orientation, and the other with the rectangle rotated 180 degrees, because you are not sure if it's oriented correctly).

5) Find which OCR output makes more sense.

6) Voila!

More Details:

1) Check this OpenCV tutorial to detect simple shapes (in your case the rectangle)
http://opencv-code.com/tutorials/detecting-simple-shapes-in-an-image/

2) if you use the previous tutorial, you have the four edges of the rectangle.
calculate the angle between the longest edge and the x-axis.

3) I don't know the level of your image processing skills, but rotating a rectangle's pixels and putting it in another image should be straight forward (I suggest using linear interpolation).

4) I don't have much experience in OCR, but according to many answers here on the exchange, I guess you can use an off-the-shelf OCR.


How can I implement a triangular filter in MatLab, given it's Z-transform representation? I have to implement an anti-aliasing filter for a certain processing step in MatLab. While searching in literature for inspiration, I came accross a paper in which the authors wrote the following:

We implement the lowpass filters as triangular smoothing. The Z-transform representation of an arbitrary N-point filter of this type is:
$$G(z)=\frac{-z^{-(k+1)}+2-z^{k+1}}{\alpha(1-z)(1-z^{-1})}$$
where the triangle length $N=2k+1$, and the filter scaling coefficient $\alpha=(k+1)^2$.

I was thinking about how to implement this filter, so following the approach taken in this DSP SE post, I rewrote the expression for the filter as follows:
$$G(z)=\frac{-z^{-(k+1)}+2-z^{k+1}}{\alpha(1-z)(1-z^{-1})}=\frac{Y(z)}{X(z)}$$
$$X(z)\Bigl(-z^{-(k+1)}+2-z^{k+1}\Bigr) = Y(z)\Bigl(\alpha(1-z)(1-z^{-1})\Bigr)$$
let's assume for now we use a filter of length $N=3$, which gives us $k=1$, we get:
\begin{align}
X(z)\Bigl(-z^{-2}+2-z^{2}\Bigr) & = Y(z)\Bigl(4(1-z)(1-z^{-1})\Bigr)\\
X(z)\Bigl(-z^{-2}+2-z^{2}\Bigr) & = Y(z)\Bigl(4(2-z-z^{-1})\Bigr)\\
X(z)\Bigl(-z^{-2}+2-z^{2}\Bigr) & = Y(z)\Bigl(8-4z-4z^{-1})\Bigr)\\
-x[n-2]+2x[n]-x[n+2] &= 8y[n]-4y[n+1]-4y[n-1]\\
\end{align}
which I finally rewrite as:
$$y[n+1]=\Bigl(8y[n]-4y[n-1]-2x[n]+x[n-2]+x[n+2]\Bigr)/4$$
I implemented this result in matlab as follows:
n = 1:1000;
x = sin(0.10*n) + (-0.5+rand(1,1000));
out=zeros(size(x));

for n_idx = 3:998
    out(n_idx+1) = (8*out(n_idx) - 4*out(n_idx-1) - 2*x(n_idx) + x(n_idx-2) + x(n_idx+2) )/4;
end

subplot(2,1,1); plot(n,x);
subplot(2,1,2); plot(n,out);

but unfortunately, the output looks absolutely nothing like a lowpass filtered version of the input. The authors do not go in detail about how they implement this filter, the only clue they give is this line:

we note that the denominator of $g(z)$ represents a causal integration of the data with $1/(1-z)$, followed by an anticausal integration with $1/(1-z^{-1})$.

which I can't make sense of.
My questions are therefore:
-1 Is my implementation of the filter $G(z)$ correct?
-2 Would anyone know what the authors mean in their comment on (anti)causal integration?

This MATLAB/OCTAVE code makes a causal implemtation the system $G(z)$ you mention. Note that $G(z)$ is non-causal, so you have to assume delays to make it causal. You can convert it to non-causal if you have the full data available before processing.

k = 5;
N = 2*k+1;

K = (k+1)^(-2);

b = K*[1, zeros(1,k) , -2, zeros(1,k), 1];
a = [1,-2,1];

[Hk,w]=freqz(b,a,1024*8,'whole');
figure,plot(w/pi,20*log10(abs(Hk)));title('frequency response');
ylabel('Magnitude (dB)');
xlabel('Frequency ( x \pi) ');

hn = impz(b,a,N);
figure,stem(0:N-1,hn);title('(causal) impulse response h[n]');


How can I create a jamming signal for desired frequency? I am working on GPS interferences and countermeasures. I came across jamming in GPS receivers.

So I have a  input signal with frequency of 1200 MHz

How can I create a jamming signal that can mask the input frequency. I am working in Matlab.

Matlab has in built function called barrage jammer, which i didn't understand.

Thanks

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

PS: Hope its only for study purposes :)


Issues related to entropy I have long been faced with the confusion regarding entropy and would be obliged if the following are answered in less technical jargon. Following the link Different kinds of entropy raises the following questions

Entropy- It is desired that the entropy of the system be maximized. Maximizing entropy means no symbol is better than the others or we do not know what the next symbol / outcome would be. However, the formula states a negative sign before the summation of the probability logarithms. Thus, it means we are maximizing a negative value!! Then if an original raw signal is quantized and the quantized information's entropy is calculated and found to be lesser than the original entropy would imply loss of information. So,why do we want to maximise entropy since it would mean that we are maximizing the uncertainty of the next symbol whereas we want to be certain about waht the next occurence of the symbol would be.
What are the differences between Shannon's entropy, topological entropy and source entropy?
What exactly is the significane of Kolgomorov complexity or Kolgomorov entropy. How is it related to Shannon's entropy?
What information does mutual information between two vectors convey?

I will try to tackle questions 1 and 4.

  1) Entropy- It is desired that the entropy of the system be maximized. Maximizing entropy means no symbol is better than the others or we do not know what the next symbol / outcome would be. However, the formula states a negative sign before the summation of the probability logarithms. Thus, it means we are maximizing a negative value!!

No, the values of the logarithms themselves are negative, so the negative sign makes them positive.  All probabilities are a real number from 0 to 1, inclusive.  The log of 1 is zero, and the log of anything less than 1 is negative.  This may seem problematic since the log of 0 is $-\infty$, but we are really trying to maximize the expected value of these logs, so when we multiply by the probability itself the entire value approaches 0, not $\infty$.  Entropy peaks when the probability is $1/2$.

  So,why do we want to maximise entropy since it would mean that we are maximizing the uncertainty of the next symbol whereas we want to be certain about waht the next occurence of the symbol would be.

No, in communicating information we absolutely do NOT want to be certain about what the next symbol will be.  If we are certain, what information is gained by receiving it?  None.  It is only through uncertainty about what the transmitter will send that we can receive any information.

  4) What information does mutual information between two vectors convey?

When there is mutual information between two vectors knowing something about one tells you something about the other.  Mathematically this equates to the following- knowledge of one vector affects the probabilities of the other vector.  If they were independent, this would not be the case.

An example of mutual information is digital walkie-talkies.  One vector is the bit stream that the first walkie-talkie sends.  The second vector is the signal that the second walkie-talkie receives.  The two are obviously related, but due to noise and unknown channel conditions the second walkie-talkie cannot know for certain what the first one sent.  It can make some really good guesses based on the signal, but it cannot know for sure.


Why do mobile phone calls produce a vibration when they are near a speaker? I remember that years ago when I had a classic mobile phone (the small one without touchscreen) when I was about to receive a call from someone and my phone was near a speaker, the speaker would start to vibrate quite hard before 2-3 seconds of receiving the call.

I don't have a way to experiment that now since I don't own a small mobile phone anymore and with modern smartphones it doesn't work. Is there any particular reason why?

How is the phenomenon produced after all? I also tried to google about it, but I had no idea how it is called.

Honestly I have no idea where I should've put this question, but searching trough all SE sites I ended up here, but if this isn't the right place please tell me.

A little background...
How a speaker works: Speakers are controlled by driving small diaphrams with a magnet and that magnet is controlled by wires that take isolated 'wired' signals and convert them into sound.

How your cell phone works: Your cell phone has antennas in it.  They work by sending and receiving transmissions of electromagnetic waves. These waves are wireless signals that are picked up by antennas and then processed later.  Your phone is transmitting/receiving those electromagnetic signals.

What is happening between the speaker and the phone: The wires in your speakers are inappropriately acting as antennas and are picking up your phones transmissions.   Note, that it only is happening when the phone is transmitting because it requires significant power to induce a current in your speaker's wires.  Also some consideration must be taken to understand the frequencies that control your speakers but that is more detail than necessary to understand the concept of what is happening.

Read this reference here if you'd like to hear more.


How does direction finding work when using GMS mobile station and additional handset? I am trying to understand the concept of direction finder (DF), for mobile GSM devices.

I have found the following description:

http://www.pki-electronic.com/products/interception-and-monitoring-systems/gsm-direction-finder/

It seems to describe the following configuration:

<IMSI catacher>-------- <additional handset (mobile)>
|
|-------- <target mobile device>

According to web site description, it is composed of IMSI catcher (probably), i.e. active base station, which force the target mobile to transmit, and probably the attacker base station (SDR radio) can detect the exact direction/signal strength of the attacked device.

I am not sure how it works, but this is what I think:

The base station (which is also a receiver) force the target mobile to keep transmission (by sending silent sms which force the mobile to transmit ack).
The additional handset is actually a second receiver.
So we get 2 receivers (handset and base stations) , each can find the signal strength.
We can draw a map, with 2 circles, and the merge points will point us to the correct direction.

Do you think this is how it works ? any comment is appreciated.

EDIT:
Can someone please give a description how it is done and what is the purpose of the handset ? Is it a reciever of just an antenna ?

Thanks.

The fact that this thing only has a single antenna port and a loudness knob:

This will just look for a transmitter in the band of interest, and you'll have to turn around a directive antenna until the tone gets loud.

That's absolutely low-tech, it can't discern between GSM phones logged into the same network (and also active), and the whole setup will only work if the straightest path to the phone is also the one that delivers the most energy – which is not the case for e.g. a lot of the indoor and urban mobile channel scenarios that we use nowadays.


$\mathcal{Z}$-transform of $\frac{1}{n^2}$ This is a Question asked in  IISC ( Indian Institute of Science,Bangalore,India) interview for MS admission.

What is the $\mathcal{Z}$-transform of  $\dfrac 1{n^2}$ ?

The problem is not sufficiently specified, because the range of admissible values of $n$ is missing. Here I make the assumption that we consider $n>0$. With this assumption we have

$$X(z)=\sum_{n=1}^{\infty}x[n]z^{-n}=\sum_{n=1}^{\infty}\frac{z^{-n}}{n^2}\tag{1}$$

And that's the point where we might get stuck, if we didn't have a list of mathematical series, or if we didn't know about the polylogarithm, which is defined by

$$\text{Li}_{s}(z)=\sum_{n=1}^{\infty}\frac{z^n}{n^s},\qquad |z|<1\tag{2}$$

where $s$ is an arbitrary complex number.

In your case, $s=2$ and the corresponding function is called the dilogarithm or Spence's function.

Comparing $(1)$ and $(2)$ we get for the $\mathcal{Z}$-transform of $1/n^2$ for $n>0$

$$X(z)=\text{Li}_2\left(\frac{1}{z}\right),\qquad |z|>1\tag{3}$$

Another way to arrive at the solution is to use the differentiation property of the $\mathcal{Z}$-transform:

$$nx[n]\Longleftrightarrow -z\frac{dX(z)}{dz}\tag{4}$$

Applying $(4)$ twice will give you the result. You need the correspondence

$$u[n-1]\Longleftrightarrow \frac{1}{z-1},\qquad |z|>1\tag{5}$$

In a first step you'll arrive at the transform of $1/n$, $n>0$, and in a second step you'll arrive at the transform of $1/n^2$, $n>0$.

In this case you will be using the integral representation of the dilogarithm:

$$\text{Li}_2(z)=-\int_0^z\frac{\ln(1-u)}{u}du\tag{6}$$


What Approximation Techniques Exist for Computing the Square Root? I have very limited resources as I'm working with a microcontroller.  Is there a taylor-series expansion, common lookup table, or recursive approach?

I'd prefer to do something without using math.h's sqrt()

http://www.cplusplus.com/reference/cmath/sqrt/

if you want a cheap and dirty optimized power-series expansion (the coefficients for Taylor series converge slowly) for sqrt() and a bunch of other trancendentals, i have some code from long ago.  i used to sell this code, but no one has paid me for it for nearly a decade.  so i think i'll release it for public consumption.  this particular file was for an application where the processor had floating point (IEEE-754 single precision) and they had a C compiler and dev system, but they did not have (or they didn't want to link in) the stdlib that would have had the standard math functions.  they did not need perfect precision, but they wanted things to be fast.  you can pretty easily reverse engineer the code to see what the power series coefficients are and write your own code.  this code assumes IEEE-754 and masked off the bits for mantissa and exponent.

it appears that the "code markup" that SE has is unfriendly with angle characters (you know ">" or "<"), so you'll probably have to hit "edit" to see all of it.

//
//    FILE: __functions.h
//
//    fast and approximate transcendental functions
//
//    copyright (c) 2004  Robert Bristow-Johnson
//
//    rbj@audioimagination.com
//

#ifndef __FUNCTIONS_H
#define __FUNCTIONS_H

#define TINY 1.0e-8
#define HUGE 1.0e8

#define PI              (3.1415926535897932384626433832795028841972)        /* pi */
#define ONE_OVER_PI     (0.3183098861837906661338147750939)
#define TWOPI           (6.2831853071795864769252867665590057683943)        /* 2*pi */
#define ONE_OVER_TWOPI  (0.15915494309189535682609381638)
#define PI_2            (1.5707963267948966192313216916397514420986)        /* pi/2 */
#define TWO_OVER_PI     (0.636619772367581332267629550188)
#define LN2             (0.6931471805599453094172321214581765680755)        /* ln(2) */
#define ONE_OVER_LN2    (1.44269504088896333066907387547)
#define LN10            (2.3025850929940456840179914546843642076011)        /* ln(10) */
#define ONE_OVER_LN10   (0.43429448190325177635683940025)
#define ROOT2           (1.4142135623730950488016887242096980785697)        /* sqrt(2) */
#define ONE_OVER_ROOT2  (0.707106781186547438494264988549)

#define DB_LOG2_ENERGY          (3.01029995663981154631945610163)           /* dB = DB_LOG2_ENERGY*__log2(energy) */
#define DB_LOG2_AMPL            (6.02059991327962309263891220326)           /* dB = DB_LOG2_AMPL*__log2(amplitude) */
#define ONE_OVER_DB_LOG2_AMPL   (0.16609640474436811218256075335)           /* amplitude = __exp2(ONE_OVER_DB_LOG2_AMPL*dB) */

#define LONG_OFFSET     4096L
#define FLOAT_OFFSET    4096.0

float   __sqrt(float x);

float   __log2(float x);
float   __exp2(float x);

float   __log(float x);
float   __exp(float x);

float   __pow(float x, float y);

float   __sin_pi(float x);
float   __cos_pi(float x);

float   __sin(float x);
float   __cos(float x);
float   __tan(float x);

float   __atan(float x);
float   __asin(float x);
float   __acos(float x);

float   __arg(float Imag, float Real);

float   __poly(float *a, int order, float x);
float   __map(float *f, float scaler, float x);
float   __discreteMap(float *f, float scaler, float x);

unsigned long __random();

#endif

//
//    FILE: __functions.c
//
//    fast and approximate transcendental functions
//
//    copyright (c) 2004  Robert Bristow-Johnson
//
//    rbj@audioimagination.com
//

#define STD_MATH_LIB 0

#include "__functions.h"

#if STD_MATH_LIB
#include "math.h"   // angle brackets don't work with SE markup
#endif

float   __sqrt(register float x)
    {
#if STD_MATH_LIB
    return (float) sqrt((double)x);
#else
    if (x > 5.877471754e-39)
        {
        register float accumulator, xPower;
        register long intPart;
        register union {float f; long i;} xBits;

        xBits.f = x;

        intPart = ((xBits.i)>>23);                  /* get biased exponent */
        intPart -= 127;                             /* unbias it */

        x = (float)(xBits.i & 0x007FFFFF);          /* mask off exponent leaving 0x800000*(mantissa - 1) */
        x *= 1.192092895507812e-07;                 /* divide by 0x800000 */

        accumulator =  1.0 + 0.49959804148061*x;
        xPower = x*x;
        accumulator += -0.12047308243453*xPower;
        xPower *= x;
        accumulator += 0.04585425015501*xPower;
        xPower *= x;
        accumulator += -0.01076564682800*xPower;

        if (intPart & 0x00000001)
            {
            accumulator *= ROOT2;                   /* an odd input exponent means an extra sqrt(2) in the output */
            }

        xBits.i = intPart >> 1;                     /* divide exponent by 2, lose LSB */
        xBits.i += 127;                             /* rebias exponent */
        xBits.i <<= 23;                             /* move biased exponent into exponent bits */

        return accumulator * xBits.f;
        }
     else
        {
        return 0.0;
        }
#endif
    }

float   __log2(register float x)
    {
#if STD_MATH_LIB
    return (float) (ONE_OVER_LN2*log((double)x));
#else
    if (x > 5.877471754e-39)
        {
        register float accumulator, xPower;
        register long intPart;

        register union {float f; long i;} xBits;

        xBits.f = x;

        intPart = ((xBits.i)>>23);                  /* get biased exponent */
        intPart -= 127;                             /* unbias it */

        x = (float)(xBits.i & 0x007FFFFF);          /* mask off exponent leaving 0x800000*(mantissa - 1) */
        x *= 1.192092895507812e-07;                 /* divide by 0x800000 */

        accumulator = 1.44254494359510*x;
        xPower = x*x;
        accumulator += -0.71814525675041*xPower;
        xPower *= x;
        accumulator += 0.45754919692582*xPower;
        xPower *= x;
        accumulator += -0.27790534462866*xPower;
        xPower *= x;
        accumulator += 0.12179791068782*xPower;
        xPower *= x;
        accumulator += -0.02584144982967*xPower;

        return accumulator + (float)intPart;
        }
     else
        {
        return -HUGE;
        }
#endif
    }

float   __exp2(register float x)
    {
#if STD_MATH_LIB
    return (float) exp(LN2*(double)x);
#else
    if (x >= -127.0)
        {
        register float accumulator, xPower;
        register union {float f; long i;} xBits;

        xBits.i = (long)(x + FLOAT_OFFSET) - LONG_OFFSET;       /* integer part */
        x -= (float)(xBits.i);                                  /* fractional part */

        accumulator = 1.0 + 0.69303212081966*x;
        xPower = x*x;
        accumulator += 0.24137976293709*xPower;
        xPower *= x;
        accumulator += 0.05203236900844*xPower;
        xPower *= x;
        accumulator += 0.01355574723481*xPower;

        xBits.i += 127;                                         /* bias integer part */
        xBits.i <<= 23;                                         /* move biased int part into exponent bits */

        return accumulator * xBits.f;
        }
     else
        {
        return 0.0;
        }
#endif
    }

float   __log(register float x)
    {
#if STD_MATH_LIB
    return (float) log((double)x);
#else
    return LN2*__log2(x);
#endif
    }

float   __exp(register float x)
    {
#if STD_MATH_LIB
    return (float) exp((double)x);
#else
    return __exp2(ONE_OVER_LN2*x);
#endif
    }

float   __pow(float x, float y)
    {
#if STD_MATH_LIB
    return (float) pow((double)x, (double)y);
#else
    return __exp2(y*__log2(x));
#endif
    }

float   __sin_pi(register float x)
    {
#if STD_MATH_LIB
    return (float) sin(PI*(double)x);
#else
    register float accumulator, xPower, xSquared;

    register long evenIntPart = ((long)(0.5*x + 1024.5) - 1024)<<1;
    x -= (float)evenIntPart;

    xSquared = x*x;
    accumulator = 3.14159265358979*x;
    xPower = xSquared*x;
    accumulator += -5.16731953364340*xPower;
    xPower *= xSquared;
    accumulator += 2.54620566822659*xPower;
    xPower *= xSquared;
    accumulator += -0.586027023087261*xPower;
    xPower *= xSquared;
    accumulator += 0.06554823491427*xPower;

    return accumulator;
#endif
    }

float   __cos_pi(register float x)
    {
#if STD_MATH_LIB
    return (float) cos(PI*(double)x);
#else
    register float accumulator, xPower, xSquared;

    register long evenIntPart = ((long)(0.5*x + 1024.5) - 1024)<<1;
    x -= (float)evenIntPart;

    xSquared = x*x;
    accumulator = 1.57079632679490*x;                       /* series for sin(PI/2*x) */
    xPower = xSquared*x;
    accumulator += -0.64596406188166*xPower;
    xPower *= xSquared;
    accumulator += 0.07969158490912*xPower;
    xPower *= xSquared;
    accumulator += -0.00467687997706*xPower;
    xPower *= xSquared;
    accumulator += 0.00015303015470*xPower;

    return 1.0 - 2.0*accumulator*accumulator;               /* cos(w) = 1 - 2*(sin(w/2))^2 */
#endif
    }

float   __sin(register float x)
    {
#if STD_MATH_LIB
    return (float) sin((double)x);
#else
    x *= ONE_OVER_PI;
    return __sin_pi(x);
#endif
    }

float   __cos(register float x)
    {
#if STD_MATH_LIB
    return (float) cos((double)x);
#else
    x *= ONE_OVER_PI;
    return __cos_pi(x);
#endif
    }

float   __tan(register float x)
    {
#if STD_MATH_LIB
    return (float) tan((double)x);
#else
    x *= ONE_OVER_PI;
    return __sin_pi(x)/__cos_pi(x);
#endif
    }

float   __atan(register float x)
    {
#if STD_MATH_LIB
    return (float) atan((double)x);
#else
    register float accumulator, xPower, xSquared, offset;

    offset = 0.0;

    if (x < -1.0)
        {
        offset = -PI_2;
        x = -1.0/x;
        }
     else if (x > 1.0)
        {
        offset = PI_2;
        x = -1.0/x;
        }
    xSquared = x*x;
    accumulator = 1.0;
    xPower = xSquared;
    accumulator += 0.33288950512027*xPower;
    xPower *= xSquared;
    accumulator += -0.08467922817644*xPower;
    xPower *= xSquared;
    accumulator += 0.03252232640125*xPower;
    xPower *= xSquared;
    accumulator += -0.00749305860992*xPower;

    return offset + x/accumulator;
#endif
    }

float   __asin(register float x)
    {
#if STD_MATH_LIB
    return (float) asin((double)x);
#else
    return __atan(x/__sqrt(1.0 - x*x));
#endif
    }

float   __acos(register float x)
    {
#if STD_MATH_LIB
    return (float) acos((double)x);
#else
    return __atan(__sqrt(1.0 - x*x)/x);
#endif
    }

float   __arg(float Imag, float Real)
    {
#if STD_MATH_LIB
    return (float) atan2((double)Imag, (double)Real);
#else
    register float accumulator, xPower, xSquared, offset, x;

    if (Imag > 0.0)
        {
        if (Imag <= -Real)
            {
            offset = PI;
            x = Imag/Real;
            }
         else if (Imag > Real)
            {
            offset = PI_2;
            x = -Real/Imag;
            }
         else
            {
            offset = 0.0;
            x = Imag/Real;
            }
        }
     else
        {
        if (Imag >= Real)
            {
            offset = -PI;
            x = Imag/Real;
            }
         else if (Imag < -Real)
            {
            offset = -PI_2;
            x = -Real/Imag;
            }
         else
            {
            offset = 0.0;
            x = Imag/Real;
            }
        }

    xSquared = x*x;
    accumulator = 1.0;
    xPower = xSquared;
    accumulator += 0.33288950512027*xPower;
    xPower *= xSquared;
    accumulator += -0.08467922817644*xPower;
    xPower *= xSquared;
    accumulator += 0.03252232640125*xPower;
    xPower *= xSquared;
    accumulator += -0.00749305860992*xPower;

    return offset + x/accumulator;
#endif
    }

float   __poly(float *a, int order, float x)
    {
    register float accumulator = 0.0, xPower;
    register int n;

    accumulator = a[0];
    xPower = x;
    for (n=1; n<=order; n++)
        {
        accumulator += a[n]*xPower;
        xPower *= x;
        }

    return accumulator;
    }

float   __map(float *f, float scaler, float x)
    {
    register long i;

    x *= scaler;

    i = (long)(x + FLOAT_OFFSET) - LONG_OFFSET;         /* round down without floor() */

    return f[i] + (f[i+1] - f[i])*(x - (float)i);       /* linear interpolate between points */
    }

float   __discreteMap(float *f, float scaler, float x)
    {
    register long i;

    x *= scaler;

    i = (long)(x + (FLOAT_OFFSET+0.5)) - LONG_OFFSET;   /* round to nearest */

    return f[i];
    }

unsigned long __random()
    {
    static unsigned long seed0 = 0x5B7A2775, seed1 = 0x80C7169F;

    seed0 += seed1;
    seed1 += seed0;

    return seed1;
    }


The meaning of the sign in a phase difference Say I have 2 signals:

A = 1kHz Sine
B = 1kHz Sine, 30° lagging the first one.

If I plot the FFT-Phase of these 2 waveforms, and get the Phase values at 1kHz, I would get:

A = -90°
B = -120°

Questions:

Should the phase difference be (B - A) or (A - B)?
What is the conventional meaning for a negative phase difference?  Lagging or leading?

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

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

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


Natural frequencies of a transfer function I have been told in the university that the natural frequencies (also called $\textit{eigenfrequencies}$) are the poles of the transfer function, however, Matlab compute them as the modulus of the pole (information taken from here).

So the real question is, are the natural frequencies complex numbers or real numbers measured in rad/s?

The natural frequency is a real number, and Matlab computes it correctly by taking the magnitude of the (complex-valued) pole.

As an example, the denominator of the transfer function of a second-order continuous-time system is given by

$$D(s)=s^2+2\omega_n\zeta s+\omega_n^2\tag{1}$$

where $\zeta$ is the damping ratio, and $\omega_n$ is the natural frequency. The poles are given by the roots of $(1)$. Assuming an underdamped system ($\zeta<1)$, the poles are given by

$$s_{\infty}=-\omega_n\left(\zeta\pm j\sqrt{1-\zeta^2}\right)\tag{2}$$

from which we get

$$|s_{\infty}|^2=\omega_n^2\left[\zeta^2+(1-\zeta^2)\right]=\omega_n^2\tag{3}$$

So the natural frequency $\omega_n$ is given by the magnitude of the pole $s_{\infty}$.

Note that the natural frequency $\omega_n$ is not the actual oscillation frequency unless there is no damping (i.e., $\zeta=0$).


Proximal Gradient Method (PGM) for a Function Model with More than 2 Functions (Sum of Functions) I am currently working in signal reconstruction. I am trying to develop an algorithm where the user can plug any constraint to the main objective function (let's say chi2, least squares). I was trying to understand the proximal gradient algorithm where you have a sum of functions $ F \left( x \right) = f \left( x \right) + g \left( x \right) $.
In the case of FISTA, it is the algorithm where $ g \left( x \right) $ is the $ {L}_{1} $ Norm of the data. However, what I am trying to do right now is to use $ g \left( x \right) $ as a sum of the $ {L}_{1} $ Norm and the Total Variation (TV) function of the data. In other words:

$$ g \left( x \right) = \lambda_1{\left\| x \right\|}_{1} + \lambda_2 \operatorname{TV}(x) $$

The question is, how can apply the Proximal function in this case?
Can I do it in iterative method where I apply the Proximal of the $ {L}_{1} $ Norm (Soft Threshodling) and then the Proximal of the Total Variation Norm?

Indeed the model for the Proximal Gradient Method (Also see Proximal Gradient Methods for Learning) is in the form of:

$$ F \left( x \right) = f \left( x \right) + g \left( x \right) $$

Where usually $ f \left( x \right) $ is convex smooth function and $ g \left( x \right) $ is convex non smooth function.

Yet the model is quite flexible and you may define $ g \left( x \right) $ any way you want.
So it can be that:

$$ g \left( x \right) = {g}_{1} \left( x \right) + {g}_{2} \left( x \right) $$

Where in your case $ {g}_{1} \left( x \right) = {\left\| x \right\|}_{1} $ and $ {g}_{2} \left( x \right) = \operatorname{TV} \left( x \right) $.

Now, like in any case of the Proximal Gradient Method you need to be able to solve the Prox of $ g \left( x \right) $:

$$ \operatorname{Prox}_{\lambda g \left( \cdot \right)} \left( y \right) = \arg \min_{x} \frac{1}{2} {\left\| x - y \right\|}^{2}_{2} + \lambda g \left( x \right) = \arg \min_{x} \frac{1}{2} {\left\| x - y \right\|}^{2}_{2} + \lambda \left( {\left\| x \right\|}_{1} + \operatorname{TV} \left( x \right) \right) $$

The above problem can actually be solved directly using various (Though slow) methods.
You could also employ methods based on Alternating Direction Method of Multipliers (ADMM) which will allow basically use the knowledge of the Prox operator for each of the functions.

But it turns out you can do even more.
Your question basically is:

$$ \operatorname{Prox}_{\lambda \left( {g}_{1} + {g}_{2} \right) \left( \cdot \right)} \left( y \right) \overset{?}{=} \operatorname{Prox}_{\lambda {g}_{1} \left( \cdot \right)} \left( \operatorname{Prox}_{\lambda {g}_{2} \left( \cdot \right)} \left( y \right) \right) \overset{?}{=} \operatorname{Prox}_{\lambda {g}_{2} \left( \cdot \right)} \left( \operatorname{Prox}_{\lambda {g}_{1} \left( \cdot \right)} \left( y \right) \right) $$

Well, In the specific case for $ {g}_{1} = {\left\| x \right\|}_{1} $ and $ {g}_{2} \left( x \right) = \operatorname{TV} \left( x \right) $ it was proved to be true in Pathwise Coordinate Optimization.

A generalization was made in the great work of Yaoliang Yu in On Decomposing the Proximal Map. He found the conditions where it holds in general. You may have a look at:

Splitting Algorithms, Modern Operator Theory, and Applications - On Decomposing the Proximal Map.
Microsoft Research - Oral 2 (Also on YouTube).
Yao-Liang Yu - Fast Fradient Algorithms for Structured Sparsity.

Simulation Code

In my MATLAB code at my StackExchange Signal Processing Q62024 GitHub Repository I wrote a script to solve the following problem:

$$ \arg \min_{x} \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2} + {\lambda}_{1} {\left\| x \right\|}_{1} + {\lambda}_{2} {\left\| x \right\|}_{TV} $$

I solved it using 2 methods:

Sub Gradient Method.
Proximal Gradient Method (PGM).

In the Proximal Gradient Method (PGM) I used the trick above where to solve the Prox of the TV norm I wrote a dedicated solver which users ADMM.

I compared the results to CVX and got this:

Indeed, as expected, the Prox method is much faster (This is even without the Accelerated Prox).

I also wrote a small script to verify the result in the article above.

Key Words: Prox of Sum of Functions, Decomposition of Prox, Composition of Prox.


Is there any information we can extract from the residual The residual in estimation is the difference between the measurement and the previous estimate. My question is when I plot this quantity, what is ,if any, information I can extract or infer from this quantity? No beneficial information I can see from the figure. Any suggestions?

If your estimate is based on a linear time invariant model then you could use the following measures for the quality of your estimate: cross-correlation between the estimate and residual, and the autocorrelation of the residual. Namely the cross-correlation should be close to zero at all time, otherwise there would still be some residual dynamics that could be modeled. The autocorrelation should be close to zero at all time except at zero, where it should be close to one. If this is the case then this would mean that the residual would be white noise, which can't be predicted/modeled.


Fourier transform's sine and cosine count as N grows I remember reading that in (discrete) Fourier transform for signals with even numbered N for length, the sine and cosine count is equal. Is this correct?

A bit of analysis:

N=1, there is only DC offset, which is a cosine wave of unlimited length.

N=2, Now, in addition to the DC offset term there is a other wave. But it's also a cosine wave, since trying to represent a sine with 2 points becomes impossible as only the 0 values of it can be sampled.

N=3, Sine can be represented as well as cos so there are 1 cos, 1 sine of same frequency and the DC offset cosine.

N=4, 3 cos, 1 sine, 3rd bin added with only cos wave possible.
...

So to me it seems that there are always going to be more cosine waves than sines. Did I make a mistake?

The interpretation that DC offset is a cosine could perhaps be debated though, has there ever been debate over this term? The fact that it surpasses the boundaries of N is a bit mind boggling.

There is no Cos and Sine. There are frequencies that represent the original signal starting with 0 rate (DC), one cycle per the N points, 2 cycles per the N points... and so on. Depends of the actual best frequency estimator of the signal the pahse of the signal is determined.


Fast pitch recognition I need to detect pitch (measure signal frequency) while the musicians play music, giving a warning if they are out of tune, but music happens to be a bit too fast for FFT (Fast Fourier Transform).

Below I try to give a technical description of the problem.

Musicians play music at 90-140 bpm. This means that there are 90-140 groups of notes each minute, up to 8 (more frequently, up to 4) notes in each group (60/140/8 = 0.0536 sec, 60/90/4 = 0.167 sec), that is, notes may change at the rate of 6-19 notes per second.

The music uses a logarithmic scale (see the attached image): the range between, say, 440Hz and 880Hz is divided into 12 notes, only 7 of which are used for melody. (Basically, they use only the white keys on the piano; when they want to shift the starting frequency, they use some of the black keys and don't use some white keys.)
That is, the frequency of each next note is multiplied by 2^(1/12) = 1.05946.

To make things more complicated, the A (La) frequency may vary from 438 to 446 Hz.  The string instruments in theory can be tuned, while the wind instruments depend on the air temperature and humidity, so the frequency happens to be re-negotiated by the musicians during the sound check.

Sometimes musicians and vocalists make errors in frequency, they call it "out of tune". They want a device that would inform them of such "out of tune errors". They have tuners, but the tuners require playing the same sound for about 1 sec before they start showing anything. This works for tuning, but does not work while the music is played.

Most likely, the tuner is doing FFT, and due to the formula $df = 1/T$ waits for 1 second to get the 1Hz resolution.

For A=440Hz, the difference in frequency between two notes is 440*0.05946 = 26.16 Hz, to get that frequency resolution, one has to use acquisition time of 0.038 sec, that is, at tempo=196bpm FFT is able to just distinguish two notes, at 98 bpm it is able to tell a 50% out-of-tune error provided that it starts acquisition at the very moment that the pitch changes. If we allow the pitch change in the course of an acquisition period, we get 49 bpm, which is just too slow. In addition, it is very desirable to be more precise about the frequency, say, detect a 25% or 12% out-of-tune error.

Is there a way to measure frequency (detect pitch) better than FFT, that is, with better resolution in less acquisition time? (At least 2 times better, ideally, 8-16 times better.)
In exchange, I do not need to distinguish between notes of different octaves, e.g. both 440 and 880 may be recognized as A. I do not need the linearity of FFT output, a logarithmic scale would be better. (Probably, more trade-offs are possible, just nothing else comes to my mind right now.)

Here's a really good drawing:

"Is there a way to measure frequency (detect pitch) better than FFT, that is, with better resolution in less acquisition time?"

yes there is.  or are.  there are multiple better ways to do musical pitch detection in real time that are far, far better than running an FFT.

consider :

Average Magnitude Difference Function (AMDF)

$$ Q_x[k] = \sum_n |x[n] - x[n-k]| $$

Average Squared Difference Function (ASDF)

$$ Q_x[k] = \sum_n (x[n] - x[n-k])^2 $$

Autocorrelation Function (AF)

$$ R_x[k] = \sum_n x[n] x[n-k] $$

note that i am playing fast-and-loose with the limits to the summation.

note also there there are no assumptions made about waveform shape or zero-crossings or other threshold crossings.  the only assumption is that when the lag $k$ is approximately a period (or two periods or some other integer multiple of the period length), $x[n]$ looks a lot like $x[n-k]$.  so the only assumption is that pitch is related to fundamental frequency of a periodic or nearly periodic (what i like to call "quasi-periodic") function.

my favorite is ASDF (and that is a thinly-veiled trade secret i just announced to everyone, but folks on comp.dsp knew that already).  these are all time domain, AMDF and ASDF look very similar and ASDF looks like an upside-down version of AF.  you are looking for nulls in AMDF or ASDF or peaks in AF which would correspond to potential period lengths of the quasi-periodic input.

here are a couple other tricks:

you can always correlate the most current $N$ samples against some $N$ samples delayed by $k$.  that way you are dealing with the most current data possible in the real-time application.
you need not calculate the correlation for every integer lag $k$.  in fact, since you like log-frequency, the spacing for larger $k$ might be bigger than the spacing for smaller $k$.
when a potential null (AM_F) or peak (AF) is found, you can compute the correlation for adjacent integer values of $k$.
between adjacent integer values of $k$, you can do interpolation to determine the peak location to a fractional-sample precision.  i won't tell you how.  use your imagination.
the whole trick (and this is the secret sauce where trade secrets and IVL patents apply) is to choose the correct peak or null when there are multiple candidates.  choosing the incorrect peak or null will result in an "octave error".  i'm not gonna tell you how to do that.  use your imagination.

send me an email and we can discuss contracting terms if you want me to design you a kick-ass pitch detector.  much better than YIN, which, in my opinion, works like shit.


Cant find Doppler shift with LFM signal Im stuck on this one. I hope someone helps me!

Problem: I need to find Dopler shift with array of reflected chirps.

I have bursts of chirps which then i sampled, and calculated fft.
Im suposed to calculate Doppler shift using just spectrum. I see it's shifted but have no idea how to calculate the shift. I tried with boxcars in hope it will  be easier and i will notice something, but it's the same thing. Also i find some really useful article and also stuck understanding it.

  "Unless the target is moving at an extremely high speed relative to
  the speed of light, the Doppler shift will be small and very difficult
  to detect from one pulse. The solution to this problem is to transmit
  a burst waveform containing repeated pulses. Initially, for
  simplicity, boxcars were used as the pulse so that some insight into
  how the parameters of the pulse affect the Doppler shift might be
  gained."
  source

The resulting spectrum appears to be your Doppler shift (possibly as a ratio to half the sampling rate, depending on what the horizontal axis on your plot represents). Without seeing all the processing that was done to generate the "Filtered Chirp Sequence" I cannot be certain but I assume it is a complex result from the subsequent plots, which would rotate in phase at the rate of the Doppler offset. To confirm this, run the test or simulation with known 0 Doppler or Clock offset and confirm that the resulting phase is constant from sample to sample.


MATLAB produces too short output for frequency modulation I have a voice data contained in a vector. Vector length is 118796 and sampling frequency is 8000 Hz. So it's around 15 seconds of voice data.

When I doing frequency modulation with fmmod (Carrier 138 MHz, out sampling frequency 5*Carrier, freq. deviation 12.5 KHz) I getting the vector of the same size.

[voice fs] = audioread('female.wav');

Fc = 138 * 10^6;
Fs = 5*Fc;
frequency_deviation = 12.5 * 10^3;
modulated = fmmod(voice, Fc, Fs, frequency_deviation);

I thing output vector should be significantly longer because it represents much higher frequency but have the same "real time" duration of 15 seconds.
Can anyone explain why matlab produce such output and how to get correct realtime output vector? I only assume to make input signal interpolation to make input vector as long as output vector should be.

The function fmmod cannot know the sampling frequency of your input signal. It assumes that the input signal has the sampling frequency Fs that you specify in the function call. Check out the example in the fmmod documentation.

By the way, in a software simulation, there's usually no need to choose such a high carrier frequency.