How to discern the distribution function of mean-squared coherence values? This probably walks a fine line between Cross Validated and this site, but I think the practical application of it lies more on the signal processing side.

I have had the occasion to record some EEG data for an experiment, and determined the mean-squared coherence between each lead:

$$ C_{xy} = \frac{\lvert{G_{xy}}\rvert^2}{G_{xx}G_{yy}} $$

(I also performed an event-related coherence calculation within certain frequency bands, which gives me a power value, but doesn't buy me anything further in the analysis of significance)

What I end up with is a list of values for the pairwise interaction of the electrodes:

Site1-Site2 = 0.24
Site2-Site3 = 0.13
...
Site7-Site45 = 0.37

so the challenge is to determine, within this sample of mean-squared coherence values, whether the MSC of the Site1-Site2 interaction is significantly different from that of the Site7-Site45 interaction, given a specific set of experimental conditions (e.g., the subject is tapping their foot or something).  Of course, ultimately comparing between subjects is necessary, but that would follow from being able to compare the results of one subject, and is probably more an issue of statistics.

Knowing that the data are not normally distributed, my question is, what is an optimal method of determining whether these values are significant?

I am familiar with 3 approaches that may be used:

Take everything above 0.5 as significant.  This is one that people could probably argue about over in stats land, but even with my untrained eye, that seems like its fudging it a bit
Use surrogate data to construct a distribution -- the signal processing piece of it would be, what would be a realistic set of surrogate data to use
Use the Bootstrap method, or basically derive the distribution from the data.  I get the idea of this method, but I don't fully understand whether this is any more realistic than using the surrogate data, or whether or not this would be SOP in signal processing

Would any of these 3 methods be appropriate for this type of analysis?  Is there another method that I've missed?

I haven't searched the literature on this in a while, but I never found a consensus anyway.  Are these standard sets of tests that you would perform on signal processing data?

One approach might be an old method from "An approximation to the cumulative distribution function of the magnitude-squared coherence estimate" by Nutall and Carter (1981), who basically apply a non-linear distortion to transform the MSC to near-Gaussian random variables, then you could apply standard methods.

http://ieeexplore.ieee.org/xpl/freeabs_all.jsp?arnumber=1163657

A more recent approach used overlapping segments to create an (approximate) PDF and CDF for the MSC: "Approximation of statistical distribution of magnitude squared coherence estimated with segment overlapping" by Bortel and Sovka (2006):

http://www.sciencedirect.com/science/article/pii/S0165168406003471

They give a procedure for the comparison of two MSC estimates.

There may be more refined methods, but I'm not aware of them.


QAM-spectrum with uneven real & non-odd imaginary part Despite the signal being real, as is seen in the last line of the first picture:

, the exemplary spectrum of the signal shows neither an even real nor an odd imaginary part:

But from my theorems I know that for every real x(t) the real part of the Fourier-Transform is even and the imaginary part of it is odd?

You are right in that the Fourier Transform (FT) of every real $x(t)$ has even real part and odd imaginary part. However, your figures only show the absolute value of the bandpass spectrum, so you cant see real- and imaginary part separately. Note that the above figure (showing $|S(f)|$ is the baseband signal, which is indeed non-even and hence complex valued. But, after upconversion, the signal becomes real and if you write a small script that does the calculation of the upconversion, you will see, the FT $S_{BP}(f)$ of the bandpass signal $s_{BP}(t)$ will have mentioned properties. I have once written a walkthrough article of up-down-conversion here which might also help in understanding this process.

To clarify a bit more in math. terms:

$\underline{s}(t)=s_I(t)+js_Q(t)$ is a complex baseband signal, which consists of the real part $s_I(t)$ and the imaginary part $s_Q(t)$.
$S(f)=\mathcal{F}\{\underline{s}(t)\}$ is the Fourier Transform of the baseband signal. This $S(f)$ is arbitrary (i.e. not even or odd) since $\underline{s}(t)$ is a complex-valued baseband signal.
$s_{BP}(t)$ is the upconverted version of $\underline{s}(t)$, and is a real-valued signal.
Hence, $S_{BP}(f)=\mathcal{F}\{s_{BP}(t)\}$ has even real part and odd imaginary part.


Why does a narrower window in a Fourier transform or STFT give poor frequency resolution? In theory of FT(Fourier transform) and STFT(Short Time Fourier Transform) it is said that
"A wide window gives better frequency resolution but poor time resolution. A narrower window gives good time resolution but poor frequency resolution"

Can somebody explain why a narrower window gives poor frequency resolution?

The other answers are good, but I thought that I would try to give a more intuitive/visual answer since I am an intuitive/visual guy.

The picture below is the plot of two tones that are almost the same frequency.  One tone is plotted in red, and the other in blue.

I generated the picture in Matlab with the following code:

tone1 = sin(2*pi*.05 * (0:99));
tone2 = sin(2*pi*.0501 * (0:99));
plot(tone1)
hold on
plot(tone2, 'r')

As you can see from both the picture and the code, the two tones are very close in frequency.  They are starting to separate a bit towards the end of this short window, but they are still very similar, to the point where a DFT could not tell them apart.

If I create the same tones with the same frequencies, only I make the window longer (5000 samples this time instead of 100) we get a much different picture.  Obviously there will be more cycles in each tone...

... but that's not the interesting part.  We see the interesting part when we zoom in at the end of the window.

We see that the two tones are 180 degrees out of phase at the end of the window, which makes them very easy to distinguish.  So how did I know to pick 5000 samples?  The difference in the two tones is $.0002\pi \frac{radians}{sample}$, so $5000$ samples $* .0002\pi \frac{radians}{sample} = \pi$ radians.

As another answer mentioned, the resolution of the DFT is the sample frequency divided by the number of samples (i.e. the window length).  I implicitly made the sample frequency 1, so if we divide 1 by 5000 (the number of samples), we get that the resolution is .0002 Hz.  Our tones actually differ by .0001 Hz, meaning that tone2 will have energy in both tone1's bin and the bin next to it (and a little bit in all the others too, but that's a different story).


How to estimate taps required for Parks-McClellan filters I have some code based on Jake Janovetz's Parks-McClellan (Remez) filter generating code. How can I estimate the number of taps required to build a lowpass filter given requirements for pass band ripple and stop band attenuation? I already know how to convert from these requirements back to the filter error deviation.

If you're using MATLAB, the function firpmord exists to help you with that. Like some other MATLAB functions, it doesn't link to any libraries or mex files, it's simply MATLAB code that runs. The only reason I mention it is that when you open this function (open firpmord), it has a subfunction remlpord that was written by (ta-daaaaa!) J. H. McClellan himself. It's using a matrix of hardcoded numbers and references Rabiner & Gold, Theory and Appications of DSP, pp. 156-7. The method therefore must be somewhat empirical, though I won't argue one way or the other. In any case, you can study this function (it's very short) and write your own based on it. I failed to find any specific papers that address the problem though, but perhaps there are books.


Hilbert transform linearity Explain why the Hilbert transform of $f(t)=\operatorname{sinc}(at) \cos(2 \pi \nu_c t)$ is

$$\hat{f} (t) = \operatorname{sinc} (at) \sin(2 \pi \nu_c t),$$

where $0<a<\nu_c.$

Attempt:

I have already found the Hilbert transform of $\operatorname{sinc}(at)$ here. Now I want to use linearity property of the Hilbert transformation. So writing $f(t)$ as:

$$f(t) = \frac{1}{2} \left( \operatorname{sinc}(at) e^{j 2 \pi \nu_c t} + \operatorname{sinc}(at) e^{-j 2 \pi \nu_c t} \right) \tag{i}$$

Applying linearity:

$$\hat{f}(t) = \frac{1}{2} \left( \frac{\sin^2(at/2)}{at/2} e^{j 2 \pi \nu_c t} + \frac{\sin^2(at/2)}{at/2} e^{-j 2 \pi \nu_c t} \right) \tag{ii}$$

But this doesn't look right. So why does the linearity not work? How can I arrive at the correct Hilbert transform expression given above?

P. S.

I know that this result follows from Bedrosian's theorem but I am not allowed to use that for this problem. I know the Hilbert transform of $f(t)$ has to be $\operatorname{sinc} (at) \sin(2 \pi \nu_c t)$ since I've already plotted its envelope and it looks correct:

Transform the signal into frequency domain

$$F(f) = \frac{1}{\left|a\right|} \left[~\frac{\Pi(f-\nu_{c})}{2} + \frac{\Pi(f+\nu_{c})|}{2} ~\right]$$

As you know the fourier domain representation of the hilbert operator is

$$H(f) = -j\text{sgn}(f)$$

Applying this transform means you just multiply $H(f)$ to $F(f)$. That gets you

$$\hat{F}(f) = -j\text{sgn}(f) \cdot \frac{1}{\left|a\right|} \left[~\frac{\Pi(f-\nu_{c}}{2} + \frac{\Pi(f+\nu_{c})|}{2} ~\right] $$

$$ \hat{F}(f) = \frac{1}{\left|a\right|} \left[~\frac{\Pi(f-\nu_{c})}{2j} + \frac{-\Pi(f+\nu_{c})|}{2j} ~\right] ~~ \text{Remember} -j = \frac{1}{j}$$

which is same as

$$ \hat{F}(f) = \frac{\Pi(f)}{\left|a\right|} \ast  \left[~\frac{1}{2j} \delta(f-\nu) + \frac{-1}{2j} \delta(f+\nu) ~\right] $$

Inverse Fourier Transformation of this gives you

$$\hat{f}(t) = \text{sinc}(at)\sin{\left(2\pi\nu_{c}t\right)}$$

Hope this helps!


Alternatives to offline Kalman filtering Recently I got into vehicle models and filtering in general and immediately faced with the following question.
I have the recorded GPS data from car driving on a highway. However, there is a substantial level of noise present, so I decided to filter this $(x, y)$ data. Currently I am using Kalman filtering (the extended one - bicycle model of the vehicle is assumed) for this purpose, but I somehow have a feeling that Kalman is more suited for online estimation/filtering. I have also tried various other approaches, including median and Savitzky–Golay filters, but they do not take dynamic model into the account which I think is the important information to use while filtering this data.
So, my question is: Is there an alternative to Kalman filter, where model dynamics can be specified, which is more suited for this purpose?
I just once again want to emphasize that I have solely measured time sequence (fixed interval) with $x, y$ coordinates and data needs to be filtered after drive is over (offline usage).
Thank you for your help and I apologize if my question is not stated well - this is my first post:)

If you used a commercial GPS receiver to obtain your tracking data, the raw trilateration data from satellites are already processed "using Kalman filtering (more common) or least-squares (less common) estimation algorithms" (see https://insidegnss-com.exactdn.com/wp-content/uploads/2018/01/marapr13-Solutions.pdf). Pay attention also to the counterintuitive statement that the article proves in the section *No Silver Bullet": the differences between least-squares and Kalman are surprisingly minor for the generic application in GPS receivers.
With this in mind, additional post-processing is quite unlikely to improve your data set, but instead your data can deteriorate as a result of this post-processing. It seems you cannot denoise your data any better, but you can analyze your route (whether the highway route passes "canyons", sharp corners or the like) and, knowing the features of your GPS receiver taken from the manufacturer's datasheet, estimate the precision of your tracking data.
The article referenced also gives a partial answer to your highlighted question about an "alternative".


Piano pitch detection I'm working on a piano tuning program and part of it requires real-time pitch detection.  Here is the scheme I have so far which works to some degree but could probably use some refinement.

I'm capturing mono, 44.1kHz, 16-bit PCM audio in chunks of 2^14 samples.  I combine the last 4 samples into a 2^16 length buffer, apply a Hann window to the buffer and run a FFT on it.  Then, I bucketize the results of the FFT in two resolutions.  First, I bucketize into 200 buckets and then run the HPS pitch detection algorithm at this granularity.  I don't need to get an exact frequency here, I just want to get close.  Then, I bucketize into 12000 buckets which gives me 1 cent resolution from 10Hz to 10kHz.  Once I know an approximate frequency from the 200 bin HPS algorithm, I search that range of the 12000 bin case for a peak to get a more exact frequency.

This seems to work okay for the notes in the middle of the keyboard.  What happens with the low notes is about 1.5s of mis-identification of the note as usually the 2nd or 3rd partial of the real note and then a correct identification of the note.

In all of the spectral plots I created to see what is going on, there is more width to the peaks that I would expect.  This width is visually somewhat consistent from the 200 bin to 12000 bin case.  I would have expected the peaks to be narrower in the 200 bin case.

So, signal processing is new to me so there may be things that are problems that I wouldn't think to ask about but in terms of specific questions, are the sample sizes sufficient for this task?  Is Hann the right choice of window?  Should I smooth the data as well before FFT?  How sensitive is HPS to the number of bins?  I was thinking that if I used a lot of bins then inharmonicity might not make partials overlap their fundamentals with the HPS algorithm's simple approach of dividing by 2, 3, 4, etc.

Similar to this thread:

Is there an algorithm for finding a frequency without DFT or FFT?

FFT isn't a particular efficient way of building a tuner. Better (and cheaper) methods include auto-correlation, phased locked loops and delay locked loops, etc..

One example is to use tracking of local maxima and minima to roughly hone in on the fundamental frequency and then use a local oscillator and phase locked loop to track this frequency precisely. This can track a moving fundamental during tuning quickly, continuously, and with great accuracy even if the frequency is low and if the fundamental is weak.


Converting a non-stationary random process into a WSS process by adding a random phase Here is an example where this method has been implemented.

We were trying to calculate the spectrum of a transmitted signal(Random signal/weighted pulse)

The auto correlation function of the pulse which is not WSS because it is dependent on t:

We associated a random variable(which is the delay) with an uniform PDF with the signal:

PDF of the delay/the variable is:

I was unable to grasp the concept behind it, probably my basics are a bit rusty.

By simply just associating a random variable (with an uniform PDF), how can we just make any random process a wide sense stationary process?
Whats the concept behind it?

It is not true that you can convert any non-stationary random process into a WSS random process by just adding a random phase. What is true is that you can make a (wide-sense) cyclostationary process WSS by adding a random phase. A common example is the PAM process mentioned in your question:

$$x(t)=\sum_ka_kp(t-kT)\tag{1}$$

The discrete values $a_k$ in $(1)$ are assumed to be random. The process $x(t)$ can be shown to be cyclostationary. The modified process

$$\tilde{x}(t)=\sum_ka_kp(t-kT+\theta)\tag{2}$$

where $\theta$ is a random variable uniformly distributed on $[0,T]$ is WSS. This random phase reflects our uncertainty about the phase of the signal, i.e., about the origin of the time axis. In this sense, it is physically motivated, and not just a nice "trick" enabling us to compute the power spectrum of $(2)$.


Implementation of wideband beamformer for planar array I would like to obtain a good reference (or references) on the implementation of a wideband beamformer for a small planar (rectangular) array comprised of 4 rows by 6 columns, for 24 elements in total.  This is a passive acoustic array with air-coupled microphones; the array does not transmit.  Each element in the array has the same spacing distance.  Given discrete input signals from each element, I would like to compute the beamformer/beamsteering output for a beam at a particular azimuth and elevation.  I need to be able to compute multiple beams to electronically steer the array. Is it possible to compute multiple beams all at the same time?

I suppose that I might need to use Fourier-domain techniques involving the FFT.  Are these techniques justified for my small array?  I've found an IEEE conference paper (http://www.curtistech.co.uk/papers/wideband.pdf or http://dx.doi.org/10.1049/ic:19980121) that may deal with this technique, but how would I apply this to a 2D planar rectangular array?  I am interested in the mathematics operating on discrete signals rather than  descriptions of the process.

Since this is an acoustic array, the elements are evenly spaced on a 4 rows by 6 columns grid.  The spacing between the elements is less than $\lambda/2$, where $\lambda$ is the minimum wavelength of the sound wave.  The wavelength range is between 3.31 cm and 331 cm.  This is a receiving array.  The target area should be in the far field of the array, but I am also wondering if it is possible to resolve signals in the near-field as well.

You need the polar response of each individual array element. That is the transfer function of each element as a function of azimuth and elevation (and potentially distance if you want to operate in the near field). For many arrays the polar function will be the same for all elements, but it doesn't hurt to check.
To calculate the polar response of the whole array you need to do the following

Create a spatial grid on which to calculate the polar response (azimuth and elevation) for the array
For each point on the grid calculate distance, azimuth and elevation of this point AS SEEN BY THE INDIVIDUAL ELEMENT. This is in essence a transformation between two different spherical coordinate systems.
Look up the single element transfer function for this direction in the single-element polar data. This may require interpolation between measured directions
Apply this transfer function to the input signal of this individual element
Make sure to properly account for the distances. That's typically done by a factor exp(-j*w*(r-r0)), where r is the distance between grid point and element and r0 and convenient normalization distance. A good choice is the distance at which the single element polar was measured.
Some results over all elements
Repeat for all points on the target polar grid


Why sinc function defined on minus side? What this relation on picture mean? Why response for a digital low-pass filter exists on "-" side of x axis? Impulse starts at 0 as i can understand and goes to infinity...
Please use phrases starting with "think of it this way..." or "it's like..."?

As it relates to linear filtering, the sinc function often comes up because it is impulse response of an ideal lowpass filter (i.e. one that has a perfectly flat passband and a transition region of zero width). In many applications, having such a filter response would be a great tool, as you could perfectly excise any out-of-band interference without disturbing a signal's passband at all.

However, there's a problem, as you noted in your question: $\text{sinc}(t)$ has infinite support! In other words, the impulse response of an ideal lowpass filter is infinitely long. That's not necessarily an issue (a lot of filters have infinitely-long impulse responses), but the big problem is that if the filter is causal (which it certainly must be in order to be realizable), then it will have infinite delay. So, if you were clever enough to build such a contraption, you would find that you would need to wait an infinite amount of time in order to see the filter's response for a given input signal. That tends to make such an idealized filter less useful in practice (even if you could build it).

With that said, while an ideal filter may not be practical, an approximation to an ideal filter can be. One simple way of designing filters, the window method, relies upon taking some filter prototype (an ideal lowpass filter, for instance) and multiplying its impulse response by a finite-length window function. The specific window function used will affect the properties of the resulting filter, but the typical goal is to approximate the prototype filter's response with something that is actually realizable.  As a result, a windowed ideal lowpass filter will still have an impulse response that looks $\text{sinc}$-like over some finite-length interval. Since a causal filter's impulse response starts at $t=0$, the main lobe of the impulse response won't occur until some later time; this is what manifests itself as filter delay.

The plot that you showed is for $\text{sinc}(t)$, which is the ideal, noncausal ideal lowpass filter's impulse response. Translating that to an implementable approximation requires truncating its response to a finite length (i.e. through windowing it) and shifting its response along the time axis so that the filter is causal (after shifting, its response would begin at $t=0$ as you expect).


Property for time invariance of system in terms of its impulse response From the textbook Signals and Systems by Oppenheimer ...

If we define $h_k[n]$ to be the impulse response to $\delta[n-k]$, then for the system to be time-invariant, we must have
$$h_k[n] = h_0[n-k]$$

I don't understand how the last equation follows. By definition of time-invariance, if the input is shifted by $k$ units, the output is shifted by $k$ units. From the given definition of $h_k[n]$, $h_0[n]$ is the impulse response to $\delta[n]$. Shouldn't the time-invariance property read
$$h_k[n] = h_0[n]$$
EDIT:
The text goes on to add:

Since $\delta[n-k]$ is the time-shifted version of $\delta[n]$, $h_k[n]$ is the time-shifted version of $h_0[n]$

But saying it is the time-shifted version does not mean equivalance. What gives?

Time-invariance means that if $y[n]$ is the response to $x[n]$, then $y[n-k]$ is the response to $x[n-k]$. So if $h_0[n]$ is the response to the impulse $\delta[n]$, then - if the system is time-invariant - its response to $\delta[n-k]$ must be $h_0[n-k]$.

If $h_k[n]=h_0[n]$ were true, then the response to any input $\delta[n-k]$ would be $h_0[n]$, regardless of $k$, which basically means that the system output would be independent of the input signal.


How to compute generator matrix from a parity check matrix? I have a parity matrix ("H") that is not in canonical form (the identity matrix is not on the right side).

I'm trying to programatically calculate the generator matrix ("G") from it.

The Wikipedia entry on Hamming codes talks about the relationship between parity check matrixes and generator matrixes:

http://en.wikipedia.org/wiki/Hamming_code

It says that H*transpose(G)=0

I thought I could figure out G by taking the Nullspace of H.  However, I end up with many fractional numbers and don't know how to use that.  The examples I've seen use gauss-jordan elimination to put the matrix into row-echelon form.  However, shouldn't I be able to do it numerically using svd or something like that?

Here's some code where I multiply a message by a generator matrix (the example is taken from Wikipedia).  I would like to get the same message when I multiply by the generator that I have calculated.

import numpy as np
import scipy.linalg

def nullspace(A, atol=1e-13, rtol=0):
    A = np.atleast_2d(A)
    u, s, vh = np.linalg.svd(A)
    tol = max(atol, rtol * s[0])
    nnz = (s >= tol).sum()
    ns = vh[nnz:].conj().T
    return ns

H = np.mat( [[1,1,1,1,0,0],
             [0,0,1,1,0,1],
             [1,0,0,1,1,0]] )

G = np.mat( [[1,0,0,1,0,1],
             [0,1,0,1,1,1],
             [0,0,1,1,1,0]] )

M = np.mat( [1,0,1] )

print "Message * generator=", M*G

GT2 = nullspace(H)
G2 = GT2.T

print "Message * calculated generator=", M*G2

With forward-error-correcting coding, one is working in a finite field, typically the field of two elements denoted by GF$(2)$ or $\mathbb F_2$. So, there are no fractional numbers
and no fancy methods such as singular value decomposition: you use bit-by-bit XOR additions of the rows of $H$ and Gauss-Jordan elimination to reduce $H$ to row-echelon form $[P_{(n-k)\times k} \mid I_{(n-k)\times(n-k)}]$.

Then, set $G = [I_{k\times k} \mid (P^T)_{k\times(n-k)}]$ and you are done. (For
nonbinary fields, use $[I \mid -P^T]$).  Note
that all arithmetic in the verification $HG^T = 0$ is also finite field arithmetic
with $1\cdot 1 =1$ and $1+1 = 0 = 1\oplus 1$ for the case of GF$(2)$.


What are the most suitable methods for correcting sudden DC offsets or baseline offset? and Why? See the image for the use-case or example The signal values should not get affected by the correction. Most methods I found were related to baseline drift or exponential baseline shift/decay. I am unable to find a method suitable for my use case.
In the image it is visible that after sudden dc offset the baseline got shifted. After correction it should stay same.

I'd do this nonlinearily, in two steps:

Detect the jump
Correct the jump

Detection sounds rather easy: a short (linear-phase) high-pass filtered version of your signal is compared to a threshold. From that, we know exactly the position of your jump; it's where the threshold was crossed, minus the group delay of the filter (which happens to be half the filter length).
Take the value of the original signal right after the jump. Shift the right-hand side of the original signal such that it ends up at the same amplitude as the sample before the jump. (A more elegant way would be definining a signal model, and then modeling the most likely signal progression, but you don't seem to have an overly useful signal, as you don't mention it.)


Dealing with bias in an FFT for frequency correction I'm working on using an FFT for frequency correction of a generic modulation signal. My algorithm basically attempts to balance the energy in the positive and negative frequencies of the FFT to be equal. One thing I'm wondering is if I have an inherent bias. Noting the following ordering returned by the bins of the FFT with N even:

After computing the FFT, my algorithm splits up the positive side into the bins $X_{up} = X[0] +...+ X[N/2 - 1]$ and the negative side into the bins $X_{low} = X[N/2] +...+ X[N-1]$ and then form an error signal as $e = X_{up} - X_{low}$
However, I realized this might be resulting in some bias because my positive frequencies are actually including DC, while my negative frequencies are including the nyquist component. So I'm wondering if the correct solution here is to drop both the DC and nyquist bins (I could add them to both the upper and lower terms but since I substract them they'd just cancel anyways) and instead do the following:
$X_{up} = X[1] + ... + X[N/2 - 1]$ and $X_{low} = X[N/2 + 1] +...+ x[N-1]$. This results in an equal number of bins, N/2-1, for both the lower and upper terms and I assume will be unbiased?

Your approach is sensible: if $N$ is even, you can discard the DC and Nyquist bins, which results in $N/2-1$ values at positive and negative frequencies. This will result in a fair comparison.
If $N$ is odd, there is no Nyquist bin, and you just discard the DC bin, resulting in $(N-1)/2$ values at positive and negative frequencies.


How to prove $cos(t) + cos(\pi t)$ is non periodic function? Also can I represent this signal using fourier series? I would just want to prove $\cos(t) + \cos(\pi t)$ is non periodic.

I don't know where to start it. Also I know that individually these signals ie $\cos(t)$ and $\cos(\pi t)$ are periodic with frequencies $1/(2 \pi)$ Hz and $1/2$ Hz respectively. But since the superpositon of these signals is not periodic should we represent them using Fourier transform only? What is the harm in drawing the frequency components at $1/(2 \pi)$ Hz and $1/2$ Hz using fourier series?

$cos(t)$ and $cos(\pi t)$ are both periodic with periods $2\pi$ and $2$ respectively. To find the period of the sum we need to find an integers $n,m$ such that $\pi n = 2 m$ or $\pi / 2 = m /n $, which is not possible since $\pi$ is irrational.

I think the Fourier series uses specific frequencies to form the signal. You can't necessarily put a frequency component at irrational frequencies.

Fourier series work fine with irrational frequencies too, but the frequencies are all harmonics, that is, integer multiples, of the fundamental frequency.


Why can the number of zeros be no more than the number of delay elements in a signal flow network? Let $N$  be the signal graph representation network of the system function (in rationale form) of a discrete-time LTI system.  For a network $N$  with
no loops, the impulse response is no longer than the total number of delay elements
in the network. From this, we conclude that if a network has no loops, then the system
function has only zeros (except for poles at $z = 0$), and the number of zeros can be no
more than the number of delay elements in the network.
Q: Why  can the number of zeros  be no
more than the number of delay elements in the network?

There is nothing preventing you from writing down an expression that has more zeros than poles.  In fact, a typical FIR filter may be designed as something like $$\sum_{n = -\frac{N-1}{2}}^{\frac{N-1}{2}} f(n) z^n$$
giving you a filter that's nice and symmetrical around zero delay.
The problem is, that in the physical world, $z$ is a lead operator -- it fetches samples from the future.  Since physicists stubbornly refuse to invent time machines for us, we have to make do with filters that only operate on past values of the input.
Thus, we only allow $z^{-n}$, and accept $z^0$ because we have an unspoken agreement that it means "the latest sample, with just a bit of delay".


What type of z transform is this? I am trying to understand this equation:

It comes from:

$Force = mass * acceleration$

$F(t) = m * a(t)$

$F(s) = m * (s^2 * y(s) - s*y0 - y0)$, where $y0=0$

$F(s) = m * s^2 * y(s)$

$y(s) =  F(s)/(m*s^2)$

Then I don't understand what type of substitution they are using for s. It appears they substitute:

$s =  (1-z^{-1})/T$

And then the last $z^{-1}$ is just to indicate they are delaying the whole thing one sample for another reason.

What type of s substitution is this called?

Have a look at these notes.

In particular, this slide:

shows that using a backward difference model sets up the approximation:
$$
\frac{d}{dt} \approx \frac{z - 1}{Tz} = \frac{1 - z^{-1}}{T}
$$

There are many ways to approximate the continuous-time derivative with discrete-time difference operations. That's all that substitution is doing: choosing one way of making the approximation.


Matlab: How do I implement the Expectation operator for time series I have the output of a time series  y(t) which is a vector  for t = 1,2,...m denoting the number of data points.

Q1: If E[y(t)] = mean of a process y. This can be invoked using the command mean(y) or
writing the statement

m = size(y,1);
Expectation_y = sum(y)/m;

Then, how
do I  find the Expectation of the following expression as

 [E[y(t-1)*y(t-1)']]^-1

These kinds of related Expressions are found in the
Slide2 mentions the Autocorrelation matrix. In the formula, there is the Expectation operator. So how do I implement the expectation of the product of the lagged random variable with itself and other such expressions without using the inbuilt commands?

My problem is how do I implement them?

Q2: My confusion is what is the different between $E[\sum_{t = 1}^ N y(t-1)*y(t-1)^T]$ and expression  $E[y(t-1)*{y(t-1)}^T]$ which looks like the Autocorrelation of the lag and again how to implement them.

Thank you!

I'm not sure I understand the question. If I'm missing something, leave a comment and I'll try to edit/expand/clarify as needed.

First of all, "mean" and "expected value" are related concepts, but are not the same thing. You find the mean of a set of samples, or observations. If you sampled/observed properly, the mean will be close to the expected value, which is the theoretical, "true mean" of the underlying random process.

When you sample a random process, you end up with a set of random numbers. If you sample the same process again, you end up with a different set of numbers. So, you can think of each sample as a random variable (a black box that produces random numbers); the random process is then a sequence of random variables.

These random variables may be correlated; that means, intuitively, that if one of them has a certain value, you can use it to predict the values of some of the other samples. For example, you may find that if one sample is positive, its neighboring samples tend to be positive too.

In your question you say that you have one set of samples, y(t). In general, you don't have enough information to generate the correlation matrix. However, you could do it like this:

M = length(y);
C = zeros(M,M);
for i = 1:M
    for j = 1:M
        C(i,j) = mean(y(i)*y(j));  % Note that the mean is useless here,
                                   % since y(i)*y(j) is scalar
    end
end

To find a closer approximation to the actual correlation matrix of the underlying process, you need to gather many observations. Let's say that you have stored the observations in the rows of a matrix Y (every row is a different set of observations). Then you could do:

M = size(Y,2);
C = zeros(M,M);
for i = 1:M
    for j = 1:M
        C(i,j) = mean(Y(:,i).*Y(:,j));
    end
end

To test this, let's try sampling a Gaussian process of zero mean and variance 1. Every time I look at the process, I'll gather three samples, and I'll look at the process 1000 times:

Y = randn(1000,3);

If I run the code above with matrix Y, I obtain this correlation matrix:

0.9598    0.0590    0.0238
0.0590    0.9675   -0.0286
0.0238   -0.0286    1.0055

which is expected (the "true" correlation matrix is the identity).

Regarding calculating a time shift in Matlab, it's helpful to arrange vectors in a table as follows (I've assumed $t=1,2,\ldots$ and $y(t)=10t$):

t    | y(t)  |  y(t-1)  | y(t-2)
     |       |          |
1    |  10   |  --      | --
2    |  20   |  10      | --
3    |  30   |  20      | 10
4    |  40   |  30      | 20
...

It's clear from the table that, to find yn=y(t-n), you cut n elements from the end of y, and shift the time vector:

tn = t(n:end);    % shifted time vector
yn = y(1:end-n);  % shifted signal

The following code provides an example, which doesn't assume integer time or integer delay:

t = 0:0.1:3;  % time vector
y = exp(t/5); % y(t)
n = 0.8;      % time delay
tn = t(t>n);  % delayed time vector
yn = y(1:length(tn));  % delayed signal
plot(t,y,'r',tn,yn,'g');  % plot

This produces the following plot:


Disparity in phase calculation between STFT and Hilbert transform I have a signal which includes several frequencies (EEG), and I want to calculate the phase of this signal in a particular frequency in a particular moment.

For this I have used two approaches a Short Time Fourier Transform (STFT) and a Hilbert transform:

I'm not worried about a possible constant (or close to constant) phase difference between the two, but about a non constant difference, which is the case, the two phases do not correlate very much.

So my impression is that both differ for the reasons below with a longer explanation on the method:

My signal has a frequency sampling of 2ms.

Let's call the time of interest T0. And foi our frequency of interest (10 hz for this particular case, though it changes).

For a particular signal segment (800 ms) I calculate an spectrogram and find the closest peak in the spectrogram to the foi, and then I calculate the a small bandwidth around it let's call it Hilbert Bandwidth (~+- 1.25Hz). (I do this to make sure there is some signal to filter later).

For the STFT,  and I use a 200 ms window centered over T0. Therefore 100 points (given 2 ms sampling) . Giving me a frequency resolution of 5 Hz. As I want the phase in a particular moment, I tried to use a window as small as possible. Then I calculate the of the foi based ont his transform.

For the Hilbert Transform I apply a bandpass Butterworth filter (order 2) using the Hilbert Bandwidth, and use a Hilbert transformation to calculate the phase.

This two phases (Hilbert and STFT) differ, so my reasoning here is:
- That the STFT phase for the foi has contributions of closer frequencies that I cannot separate due to the length of the window. While the Hilbert does not as I have applied a bandpass filter narrower than the resolutuion of the STFT.

So I have three questions:

Is filtering in a narrow band and calculating the Hilbert transform a senseful way of calculating the phase of the signal?
Does my explanation of the difference between the two methods make sense?
For a given signal,  which of the two methods would be more reliable for estimating the phase at a particular time point in a particular frequency?

Without seeing your data, it's hard to tell whether you have contributions of closer frequencies or not.  That is certainly possible.  So, my suggestion to you would be this: try out your two different methods on pure sine waves where you know that there are no closer frequencies.  In other words, generate some dummy data which is just one single frequency.  If your two methods give the exact same answer on that data, then maybe it is in fact the contributions of closer frequencies that are messing up your original calculation.  But if the two methods do not give the exact same answer, then you likely have an error in your method.


How to deal with the mirrored components of a FFT? And another question I'm trying to simulate the behavior of my FFT spectrum analyzer, I am using Mathematica, but my questions are more conceptual than technical implementation of code...I hope.

If we consider a simple cosine wave transient signal $y(t) = A \cos(2 \pi \nu_0 t)$, where $A = 1$ and $\nu_{0} = 24$ $\rm{Hz}$. This signal is sampled at some time interval $dt$ and this gives me 2048 points (2048 because of my device), I want to resolve this data to some frequency resolution, $\Delta f$, and from the Fourier limit: $t_{sample} = 1/\Delta{f}$ so my sample rate, $S_{r} = 1/dt$ is adjusted accordingly.

I then take the absolute value of the FFT of these 2048 points and reconstruct the frequency component using $\Delta f$ (in this example $\Delta f = 0.0625$ $\rm{Hz}$, this value is because of my device and the features I want to resolve) and the index of the resultant FFT, and get this:

I understand that the symmetric appearance of an FFT comes from the real and imaginary components having the same response to signals -- the absolute value of their responses will be the same. But the phase of the real and imaginary parts is shifted by $\pi / 4$. This gives us our mirrored appearance.

FIRST QUESTION: Usually I would just throw the right hand frequencies away and multiply the amplitude of what is left by $2$. Is this the correct approach or should I reverse the order of one half of these 2048 points and then add them together (intuitively I would say no, because both halves should be identical so multiplying by $2$ should suffice)? What is the best method of reconstructing the frequency component/axis?

SECOND QUESTION: This is something I just don't understand. If I increase the frequency of my transient signal $\nu_{0}$ then the two mirrored peaks move closer to each other, converging in the center. This makes sense. If I keep increasing $\nu_{0}$ of the transient arbitrarily (past $124$ $\rm{Hz}$ which corresponds to the limit of frequencies) however I still see a peak, in fact I see this oscillating behavior of the peaks just bouncing back and forth as I further increase the frequency -- Why?! Is this natural or because of the way I define my frequencies.

It depends on your spectrum analyzer.

If you want to analyze a one channel real baseband signal, the mirrored frequencies are redundant.  you should keep the first N/2+1 points.

For communications signals, there can be an I and Q (two reals as a single complex), you want to keep them all.  the frequencies will not be (in general) symmetric.

If you want to calculate the cross spectrum of two reals, there will be mirrored response, but if you want to do something like taking the inverse of the cross spectrum, you should keep both, but only display half.

It really depends on the actual application. There really isn't a right or wrong way. There are conventions and exceptions. There are niche applications. Spectrum analyzers that are commercial standalone test equipment, usually have lots of buttons, lots of modes.

for your second question, this is what spectral aliasing looks like if you don't apply  an antialiasing filter prior to the DFT.


The Bode frequency shifter Anyone know if the Bode frequency shifter has been studied somewhere. I just used an audio plug-in called EchoBode from SonicCharge and this is claimed to use a frequency (pitch) shifting algorithm inspired by the Bode frequency shifter.

The pitch shifting quality on this is pretty darn good I think, thus why I'm asking about it.

Here it's claimed to involve a ring modulator and a Hilbert transform.

Frequency shifting is linear and thus produces inharmonicity, but for small (max. 1-2 semitones) the quality is possibly better than with "granular" pitch shifting approaches. The inharmonicity is also not necessarily distortion or "ring-modulated", but may be even musical for some sounds such as percussive instruments.

They may be referring to patent US3800088 by Harald Bode.

I have a bunch of images from 15 years ago to explain a way to do frequency shifting so that it has the best chance of sounding good. I would call it single-sideband frequency shifting. Here the range of possible frequencies in the signal are drawn as two arrows each spanning half the circumference of the z-plane unit circle:

T is integer and represents time. The frequency shift $-\frac{\pi}{4}$ is an example. Modulation here is by a complex exponential (the signal is multiplied by the complex exponential).

The resulting signal is complex, but we want a real signal, which we get by throwing away the imaginary part. This creates negative frequencies (green) that are mirror images of the positive frequencies.

That's pretty ugly as some frequencies got shifted so that they alias and thus have a positive frequency shift. Also the spectrum of what we want got mixed with wrong-way-going frequencies (reversed spectrum). We can fix these problems using the following procedure:

We filter away negative frequencies with a filter that has a do-nothing real part and a Hilbert transformer as the imaginary part. Then we modulate with a complex exponential and use another positive frequency pass filter to remove anything that might alias. Finally, discarding the imaginary component creates a nice-looking signal.

I excluded details about gain, but it suffices to know that if at the end the gain is wrong, the output needs to be multiplied by some factor $2^N$, $N$ integer.


CT corrupted by Poisson noise, weighted least squares, how to calculate in detail I've trained a convolutional neural network to recognize ellipses. I took random generated ellipses as ground truth images and applied radon and inverse radon transform for the noisy images.

Now I want to add Poisson noise to the sinogram before applying the inverse radon transform.
For the reconstruction, the paper I'm following suggests to replace the data-fidelity term in

${min}_x \frac{1}{2} || Hx-y||^2_2$

with weighted least squares. (Here, x is the input image and y is what we are measuring).

It goes on with: "Let $\mu$ represent the distribution of linear attenuation coefficient of an object and $[A \mu]_m$ represents their line integral. The mth CT measurement, $y_m$ is a Poisson random variable with parameters

$ p_m \sim Poisson( b_m e^{-[A \mu]_m} +r_m)$

$ y_m = -log(\frac{p_m}{b_m})$

where $b_m$ is the blank scan factor and $r_m$ is the readout noise."

The paper then suggests to use these $p_m$ to build a diagonal weight matrix $W$ to get

${min}_x \frac{1}{2} || W^{1/2}Hx-W^{1/2}y||^2_2$.

But what exactly is $[A \mu]_m$? Where do I get $A$ and $\mu$?
In a linked paper, it says

$[A \mu]_m=\sum_{j=1}^J a_{mj} \mu_j, m=1, \dots , M$, where $A=\{a_{ij}\}$ is the system matrix and M and J is the the number of projections and pixels.

But again, what exactly is the system Matrix $A$, where can I find it?

Here are the papers I'm referring to: number 1 (page 12 (1451), part B) and number 2 (page 3 (259), part B).

Any help is highly appreciated, thank you very much.

Edit:

Thanks to @A_A, I now have an idea of how $A$ looks like. But I still don't really know how to calcalute $p_m$. Would it be possible to get $p_m$ from $y_m=−log(\frac{p_m}{b_m})$, that means $p_m=b_m e^{−y_m}$? But how exactly would I determine $y_m$; let's say my image $x$ is from $\mathbb{R}^{512×512}$ and performing a radon tranformation gives me a sinogram $Hx \in \mathbb{R}^{512×32}$. Then $W \in \mathbb{R}^{512×512}$ and there are $512$ values of $p_m$. How would I go on from there? And what should I choose for the blank scan factor?

Thank you very much.

The first paper is more clear than the second.

"Let $\mu$ represent the distribution of linear attenuation coefficient of an object and $[A \mu]_m$ represents their line integral. "

The signal obtained from Computed Tomography (CT) depends on how much was the X-Ray beam attenuated as it travels through the human body (there are other modalities such as the Positron Emission Tomography (PET) where the signal depends on how much emission you measure. Same reconstruction technique, different signal model).
You don't really know what is inside some $A$. It represents the body. What you do know is that it comes in the form of a matrix of sites which is called $A$. Now that you have the $A$, then you consider the total attenuation along the path of a ray through this $A$. That is the "line integral" $[A \mu]_m$ .

The $m^{th}$ CT measurement, $y_m$ is a Poisson random variable with parameters
$ p_m \sim Poisson( b_m e^{-[A \mu]_m} +r_m)$
$ y_m = -log(\frac{p_m}{b_m})$
where $b_m$ is the blank scan factor and $r_m$ is the readout noise."

This $[A \mu]_m$ is inaccessible to you. What is accessible to you is $y_m$, which is what you measure, it contains noise and is expressed here as the ratio of a measurement through the body divided by a measurement in air ("blank scan factor").
Notice here, $y_m$ is the line integral. One number. One ray through $A$. That is, one pixel from one line of your Radon transform. One measurement. NOT the whole line of the Radon transform. To look along the line of the Radon you need to do it in an iteration (i.e, expressions where $\sum_{m=1}^M y_m \times ...$).
Hope this helps.
EDIT:

I still don't really know how to calcalute $p_m$.

$p_m$'s expression is given in the paper as coming from a Poisson distribution. It does not look very much Poisson because it is missing the rate ($k$) and so it ends up being closer to an exponentially decaying distribution. The effect this is trying to simulate is the variation of the detector's sensitivity with attenuation (or radiation energy).

Would it be possible to get $p_m$ from $y_m=−log(\frac{p_m}{b_m})$, that means $p_m=b_m e^{−y_m}$?

That means $p_m \sim Poisson( b_m e^{-[A \mu]_m} + r_m)$. $r_m$ is the readout noise.

But how exactly would I determine $y_m$;
[...]
And what should I choose for the blank scan factor?

I don't know what should you choose for $b_m, r_m$, it depends on what you are modeling (detector voltage? Hounsfield scale?). Usually, $A$ is supposed to contain Hounsfield scale values directly, in which case your blank scan factor should be that for air. But this means that the rest of your image values are adjusted according to what sort of tissue they depict. The point of this is to take into account the Poisson noise because it modulates the readings. So, perhaps you can lump things together into reasonable values obtained by measurements over phantoms.

let's say my image $x$ is from $\mathbb{R}^{512×512}$ and performing a radon tranformation gives me a sinogram $Hx \in \mathbb{R}^{512×32}$.

Just a note here, if you keep your "detector width" at 512 measurements, you will be sampling a disc inscribed within the $512 \times 512$ square. When the detector is looking at a $512 \times 512$ image at a $ \frac{\pi}{4}$ bearing, its diagonal is $\approx 724$ pixels. What I like to do is to center the square image in a black frame that is $diagonal \times diagonal$ in dimensions and use that for reconstruction. It is also easier to keep the detector geometry fixed and rotate the image. For a parallel ray detector, the line integrals are simply the sums of the columns of the image. For a fan beam, each ray deviates but you don't have to re-calculate this for every projection.

Then $W \in \mathbb{R}^{512×512}$ and there are $512$ values of $p_m$. How would I go on from there?

As far as I can tell, this $W$ term is $w_m = \frac{(p_m - r_m)^2}{p_m}$ (equation 29).
But, the part of the discussion on page 1451 that involves $W$ is to explain how to calculate your $E$ that drives the optimisation detailed in equation 33. between what you "guess" $Hx$ (see "data fidelity term" pages 1440-1441 around equations 1 and 2) and what you measure $ym$.
So how image is affected by Poisson noise is one thing and calculating that error term during optimisation is another.
Hope this helps.


Why is there no FPU on (most) DSP chips? Why is there no FPU (Floating Point Unit) on most DSP chips?  I found
several answers such as this
one
that explain advantages of fixed point approach such as smaller power
consumption, higher speed and I can myself think of several other
advantages such as smaller heat production or price but I have never
seen any practical real life benchmarks. The questions I would like to
know answers for are:

how much is FPU and what is the price difference between CPUs with an FPU and CPUs without an FPU?
how much heat can FPU produce?
how much power would FPU consume?

An 8086 used less than 30k transistors. The 8087, which is the FPU for the 8086, is reported to use 45k transistors.  Faster FPUs can be even larger in terms of gate count.  So the cost in silicon die area of an FPU can be significant (over 2X?).  Power and thus heat is proportional on the order of the number of transistors toggling outputs at similar rates.

For real-time DSP, there is also the issue of deterministic latency. A simple integer multiply-accumulate unit has a fixed latency.  A really simple floating point adder does not have a fixed latency (or is either not fast or not simple) due to the potential need to normalize up to twice and/or handle NaNs and denorms.


Derivation of range migration algorithm Problem:
In Walter G.Carrara's book on synthetic aperture radar, the equation is presented:
$\Phi(K_X, K_R) = -K_XX_t - R_B\sqrt(K_R^2 - K_X^2) +K_RR_S $  (10.30)

And this is said to come from substituting the value for $X_a$ using the Principle of Stationary Phase:
$X_a = -R_BK_X/\sqrt(K_R^2 - K_X^2) + X_t $ (10.28)

into the phase of the fourier transformed signal:
$\Phi(K_X, K_R) = -K_RR_t - K_XX_a + K_RR_s $ (10.29)

What I have tried:
Substituting 10.28 into 10.29 as mentioned, I got:
$ -K_RR_t + K_X^2R_B/\sqrt(K_R^2 - K_X^2) - K_XX_t + K_RR_s $

The following identities may be useful:

$R_B = \sqrt((Y_a - Y_t)^2 + (Z_a - Z_t)^2)$ is minimum range to scatterer

$R_t = \sqrt((Xa - X_t)^2 + R_B^2) $


Using fft in matlab to retrieve input signal magnitude I am trying to understand how you use an FFT in matlab to measure/retrieve the magnitude of the input signal.

For my understanding I have put this script together, creating a 15Hz sine wave with a sampling frequency of 1kHz.

fftsize    = 2048;            % How many samples for FFT
fs         = 1000;            % Sampling frequency (Hz)
sim_length = fftsize/fs ;     % run sim for this time
t          = 0:1/fs:sim_length-1/fs;
x          = (0.5)*sin(2*pi*15*t); % 15 Hz component, 0.5 magnitude

Window and FFT:

%% Windowed FFT
data                  = x'; % Assign input data
window                = blackman(fftsize);
windowed_data         = window.*data;

% Perform FFT and take magnitude (of complex result)
fftd_data             = abs( fft(windowed_data, fftsize) );

The input frequency was not coherent, add in the windowing and there is going to be some spectral leakage, therefore I need to sum the power over a few bins around where I expect the fundamental to be.

% Each bin represents
bin_freq = fs/fftsize ;

% 15 Hz fundamental is centred in bin
fundamental_bin = ceil( 15/bin_freq );

% Allow spectral spreading of +-4 bins
lower_limit = fundamental_bin - 4;
upper_limit = fundamental_bin + 4;

Sum of values from fft:

sum( fftd_data(lower_limit:upper_limit) )

ans =

  512.3738

Which is not the 0.5 I was looking for. What is the best way to get the (0.5) input magnitude from the fft result.

After some digging about I have come up with the following but not quite sure why it works or if I have named the calculations appropriately.

incoherent_power_gain = sum(window.^2);
power_spectrum        = fftd_data.^2 ;
power_fftd_data       = (4 * power_spectrum) / (incoherent_power_gain * fftsize);

sum_fundamental_power = sum( power_fftd_data(lower_limit:upper_limit) )

sum_fundamental_power =

    0.2500

We converted to Power so we can take the square root to get back to a signal level.

sqrt( sum_fundamental_power )

ans =

    0.5000

Great, but is there a better way or can some one explain:

If power spectrum is the correct name for fftd_data.^2.

I am also unsure what is happening in (4 * power_spectrum) / (incoherent_power_gain * fftsize).

I just fixed 3 things in your code and added comments where I did. The reason is basicly just the way the FFT is defined/implemented in Matlab.

fs         = 1024;
fftsize    = 2048;           % Choose as a power of 2 of the sample rate

sim_length = fftsize/fs ;
t          = 0:1/fs:sim_length-1/fs;
x          = 0.5*sin(2*pi*15*t);

data                  = x';
window                = blackman(fftsize);
windowed_data         = window.*data;

fftd_data   = abs( fft(windowed_data, fftsize)/length(x) ); % Normalize the fft

bin_freq = fs/(fftsize) ;

fundamental_bin = ceil( 15/bin_freq );

lower_limit = fundamental_bin-4;
upper_limit = fundamental_bin+4;

2*sum( fftd_data(lower_limit:upper_limit) ) % Multiply by 2, because of the DFT-symmetry

ans =

    0.5000

For a more detailed explanation, check out this thread at Mathworks forum.


Simple problem regarding rectangular pulse shifting So i have the following rectangular pulse:

$$h(t)=\begin{cases}
      1 &  0<t< T \\
       0 & \text{otherwise}
   \end{cases}$$

Now, i am supposed to create $h(T-t)$, and the problem is that i do it two different ways, and get two different results, and it is obvious that this is not possible. Anyway, here is what i have.

My first attempt is step by step building of the $h(T-t)$ from the $h(t)$, firstly, i "flipped" variable $t$.
so i got
$$h(-t)=\begin{cases}
      1 &  0<-t< T \\
       0 & \text{otherwise}
   \end{cases}=\begin{cases}
      1 &  -T<t< 0 \\
       0 & \text{otherwise}
   \end{cases}$$

Now, i add T to this, and i have
$$h(T-t)=\begin{cases}
      1 &  -T<T+t< 0 \\
       0 & \text{otherwise }
   \end{cases}=\begin{cases}
      1 &  -2T<t< -T \\
       0 & \text{otherwise}
   \end{cases}$$

However, i've done this directly:

$$h(T-t)=\begin{cases}
      1 &  0<T-t< T \\
       0 & \text{otherwise}
   \end{cases}=\begin{cases}
      1 &  -T<-t< 0\\
       0 & \text{otherwise }
   \end{cases}=\begin{cases}
      1 &  0<t< T \\
       0 & \text{otherwise}
   \end{cases}$$

Now, there is something wrong here, but i don't see any mistakes in any of the two cases. Any help appreciated!

Clearly, in this example you have $h(t)=h(T-t)$. The easiest way to see this is to check for which value of $t$ the argument becomes zero. So the point $t=0$ of the original function is mapped to the point $t=T$, and because you have $-t$ in the argument, the function is inverted on the time axis.

The mistake in your first approach is the way you added the value of $T$. What you should have done is replace $-t$ by $T-t$ (which is what's happening on the left-hand side of the equation), to obtain $h(T-t)=1$ for $0<T-t<T$, which is equivalent to $0<t<T$, corresponding to the correct solution.


How to determine the sine Fourier coefficients of discrete data? The following relation gives me the measurements of interest $w$ at equally distanced locations $x_j$ in space:
$$w_j=\sum_{m=1}^{11}A_m\sin\left(\frac{mπx_j}{L}\right)$$
where $A_m$ are the Fourier coefficients of sine series and $L$ the total length is physical space. I also assume $m = 1,2,...,11$.
Now, how can I obtain the coefficients $A_m$ given the data $w_j$ ?

In accordance with the case of a continuous function $f(x)$ where
$$b_n=\frac{2}{L}\int_{0}^Lf(x)\sin\left(\frac{mπx}{L}\right)dx$$
I came with the following expression:
$$A_m=\frac{1}{N-1}\sum_{n=1}^{N}w_n\sin\left(\frac{mπx_n}{L}\right)$$
Can anyone confirm or not that solution?


Method for finding the best non-increasing approximation function of a given signal I have a discrete-time signal (with noise) and I want to determine the best non-increasing function that represents it.

Of course it would be useful if the obtained function were defined in terms of a continuous variable of time, but to me this sounds too general. I think it is not a good idea to specify a particular equation with coefficients to be determined, simply because in so doing I am imposing a too restrictive space of test functions, thus possibly causing "underfitting". For example, I could impose an exponential function like $f(t) = c \, \exp(-\lambda t)$, where $c$ and $\lambda$ are parameters to be determined, but what if the given signal has another form?

For a discrete signal I think it is better try to generate a discrete approximation as well. In this case, the "fit" would be a non-increasing discrete function that is the best representation of the given general signal.
Does anybody know a method to do that?

One could use either the raw signal or then a low-pass filter firstly. As a remark I would like to say that polynomial interpolation simply does not work because of the noise.

You will find some pointers with the correct wording. There is a lot of statistical literature on monotone or monotonic regression (sometimes called isotonic regression). A more generic term is "shape-constrained estimation".

For instance, a few references:

Constrained statistical inference: inequality, order, and shape restrictions, 2005, Mervyn J. Silvapulle and Pranab Kumar Sen
Polynomial algorithms for isotonic regression, 1997, Victor Chepoi and Daniel Cogneau and Bernard Fichet
Fitting monotonic polynomials to data, 1994, Douglas M. Hawkins

The most common version is a least square fit with piecewise constant or linear functions, see for instance Scikit-learn: Isotonic Regression

  An illustration of the isotonic regression on generated data. The
  isotonic regression finds a non-decreasing approximation of a function
  while minimizing the mean squared error on the training data. The
  benefit of such a model is that it does not assume any form for the
  target function such as linearity. For comparison a linear regression
  is also presented.

With more recent progress in optimization, there are works on fitting under other norms ($\ell_1$, $\ell_\infty$) and contraints on the lower or upper bound of the derivative, which should be either positive or negative, to impose non-increasing or non-decreasing approximations.

[EDIT:2018/12/24] Now, on the question whether the signal should be low-pass or not: if you think that a paramteric function is  too restrictive, a low-pass filter would probably too restrictive as well. Indeed, how do you know that the characteristics or the filter is sufficient? wgy should it be linear? does it cope with the nature of the signal and noise well enough?


What is Shannon's source entropy Suppose that ${X_n; Y_n}$ is a random process with a discrete alphabet, that is, taking on values in a discrete set for $n$ data length. They correspond to the input and output of a communication process. Assuming Y to be the discretized output, what is meant by Shannon's source entropy for a given discretization bin $k$ : $H_{source} = lim_{k \rightarrow \infty} \frac{1}{k} H_k$ where
 $H_k$ stands for Shannons entropy. I found this in the paper. Now, problem_Set also mentions about source entropy but the formula is very different! Its the same as Shannon's entropy. So,what exactly is Source entropy? Is it the entropy of X or Y or is it the same as Shannon's entropy

So,what exactly is Source entropy?

When the paper talks about "source entropy" all they mean is the entropy of the information source.  You can see that from the following passage in the paper- "Shannon showed that there must exist at least one encoding of the sequence generated by an information source which allows error-free transmission of the sequence when the channel capacity is larger than the source entropy."

  Now, problem_Set also mentions about source entropy but the formula is very different!

Yes, and no.  The formulas appear different, but they are not really.  The first formula is
$H_{source} = lim_{k→∞}\frac{1}{k}H_k$

$H_{source}$ is the total entropy of the source, and $H_k$ is the entropy of each possible discrete value that the random variables in $H$ can be.  Thus, it is almost a trivial statement that the total entropy is the sum of the individual value entropies, while the number of discrete values can go to $\infty$.

The problem set equation is-
$H(s) = \sum_i p_ilog_2\frac{1}{p_i}$

This is really the same equation, it just replaces the generic "$H_k$" with the implied $log_2\frac{1}{p_i}$.  The $\frac{1}{k}$ in the first equation is equivalent to the $p_i$ in the second equation- i.e. they are assuming that all of the discrete values of H are equally likely.  This is not generally true, but is true when the entropy is maximized.


Time period to Instantaneous frequency I have simple smooth time signal with different time period. I would like to transfer period over the time to frequency, and to track frequency changing over the time, I would call it Instantaneous frequency?

PF-period frequency (1/T)... any suggestion?

You can compute the instantaneous frequency of your data in Matlab as follows:

z = hilbert(x);
N = length(z);
f = angle(z(2:N).*conj(z(1:N-1)))/2/pi;

To get units of Hertz, multiply f by the sample rate.

There are other, similar estimators that do a similar thing. This matlab code can help.

Eric Jacobsen also has a couple of pages here and here.


How Random Forests can be used for regression? I have an understanding of how random forests can be used for regression. But, I am still unable to understand, how random forests can be used for regression.
I have already read the research work of Leo Brieman on random forests. But, I still fail to understand how random forests can be used for regression.

It would be of great help, if this concept can be explained using any toy example.

Thanks in advance.

The key is in choosing the split function of the decision tree. In classification, the data is split according to a certain homogeneity metric. On the other hand, for a regression tree, a regression model is fit to the dependent variable, for each of the independent variables. Then, the data is split at multiple points and the SSD value is evaluated. The best split is chosen for the dimension (variable) giving the minimum error. The node is then split at this point. The splitting is again recursive.

Finally, each leaf stores a distribution for the continuous output variable/s. In the testing phase, the responses are averaged.


What do the function blocks mean? What I don't understand is what exp(.) and ln(.) mean. For example y(t) and x(t) have exp(.) between them. What exactly does that mean? y(t)=?

Hint:
Write the relationship between $z(t)$ and $x(t)$ and then simplify it a bit.


Is mean squared error or absolute error on two spectrograms of two different wav files enough to compare how similar they are? I've been working on a genetic algorithm to figure out how to play a guitar chord just by listening to it. I have it all working now, but the system I am using to compare how 'similar' a guess is to the target seems less than robust. Here is the basic workflow of the algorithm:

I have a wav file and a corresponding spectrogram of the 'target' chord which I am trying to have the algorithm recreate
I create a population of 'chords'. I simulate how a person would place their fingers on a fretboard and then convert it to a more readable set of pitches.
Then, I run each chord through the simulation. I play what that chord would sound like through a MIDI library and, at the same time, record the sound to convert into a spectrogram. This is not ideal because then I have to wait for each chord to play in isolation, but to my knowledge, there is no way to just dream up a spectrogram, I have to just record it.
Then I assign each chord a 'fitness' value by doing either mean squared or absolute error pixel-wise between the target spectrogram and each chord's spectrogram. I also cut any silence before or after the wav file before turning them into a spectrogram and take the average across the time dimension for each spectrogram (so time is not a factor, and instead of being a plot of frequencies vs time, it is just a list of frequencies and their corresponding amplitudes)

However, I have noticed that this algorithm for finding the fitness of each organism doesn't work great. Even when I run the same chord shape through multiple times, I get a lot of fluctuation each time as to what fitness it receives. It should be relatively the same volume each time, so I can't imagine normalization will help, but I could try it. Any other ideas about how to get it more reliable? Are spectrograms not the way to go? Or maybe there is a smarter option than pixel-wise error? Thanks!

My suggestion would be to partition the spectrum according to the well tempered scale, one frequency band for each half tone. You generate a pattern of 1s and 0s for each chord you want to be able to test, with a maximum of six 1s.
When a chord is played, just do a FFT of the whole thing, normalize the energy and check each band for a threshold.
Your "fitness" value then could be directly calculated from the correct strings vs. wrong strings in some way.


Running the BM3D Denoising Algorithm in Python I am looking for a strong baseline in image denoising and therefore wanted to have the BM3D algorithm in my benchmark.

These 2 python implementations:

https://github.com/ericmjonas/pybm3d
https://github.com/liuhuang31/BM3D-Denoise

have been unmaintained for a long time (and the second one is documented in what I think is chinese). I have experimented a bit with the first one without satisfactory results (see https://github.com/ericmjonas/pybm3d/issues/11 with non satisfactory fixes like clipping).

This package : https://pypi.org/project/bm3d/#description , doesn't have a documentation or source code easily findable.

Then I found this: https://docs.opencv.org/master/de/daa/group__xphoto.html#ga2fc5a9661c1338a823fb3290673a880d, in the docs of opencv but there is no indication as to how to use it and no examples in python. I saw this SO question but it's for a C++ implementation and I don't know how it would translate in Python because I am not familiar at all with opencv.

What I resorted to was using the PyPI package, which is advertised here: http://www.cs.tut.fi/~foi/GCF-BM3D/index.html#ref_software.

I dug a bit in the source code, and found that I could perform BM3D, in the following fashion:

import bm3d

denoised_image = bm3d.bm3d(image_noisy, sigma_psd=30/255, stage_arg=bm3d.BM3DStages.HARD_THRESHOLDING)

There are also some examples in the library's source code download.

I installed bm3d using pip (pip install bm3d) and needed OpenBlas (sudo apt-get install libopenblas-dev).


Instantaneous frequency vs fourier frequency Lets consider a pure sine signal at $\nu$ that is chopped using square pulses (like a burst mode on signal generators). My understanding is that instantaneous frequency is $\nu$ when oscillations are ON and 0 when they are OFF. On the other hand fourier spectrum is constant over time and contains also other frequencies, since it is not pure sine anymore. Is this correct? which one is used when calculating some frequency dependent physical quantity?

the frequency of a sinusoid is normally thought to be undefined when there is no sinusoid (i.e. when it is turned "OFF").

now, there are at least a couple of different methods of turning it OFF.  the most obvious is to effectively multiply the waveform output by zero to disconnect the sinusoid from the output.

however if the output of your signal generator is always connected (so there is no multiplication by zero) to the look-up table (LUT) that has the waveform values in it (and let's say that those waveform values follow a sinusoidal function), then another way to stop the output is to stop advancing position in the LUT.  that would be like setting the instantaneous frequency to zero.  but, unless you're watching out for this, it could get stuck on a DC value that is non-zero.  so there are a lot of functions, besides the zero function, that have a frequency of zero.  any DC function has instantaneous frequency of zero.


Deconvolution (Linear Convolution) with an Under Determined System of Equations? If I have a measured signal $\mathbf{y}$ which is the result of the true signal $\mathbf{x}$ convolved with another function (linear and not circular convolution), I always seem to get an underdetermined system. Convolution expressed as a matrix equation:

\begin{equation}
\mathbf{y} = \mathbf{A} \mathbf{x}
\end{equation}

I am interested in finding $\mathbf{x}$, given $M=5$ observations and knowing that the convolution matrix is defined by a boxcar function of window length $n=3$. The corresponding $M\times N$ convolution matrix is defined as the following, where $N=M+n-1$.
\begin{equation}
\mathbf{A} = \begin{pmatrix}
1 &   1 &   1 &   0 &   0 &   0 &   0 \\
0 &   1 &   1 &   1 &   0 &   0 &   0 \\
0 &   0 &   1 &   1 &   1 &   0 &   0 \\
0 &   0 &   0 &   1 &   1 &   1 &   0 \\
0 &   0 &   0 &   0 &   1 &   1 &   1 \\
\end{pmatrix}
\end{equation}

In the literature, the convolution matrix is often $N\times N$ (or $M\times{}M$, however you want to define it). But it seems that with a finite number of observations, the convolution matrix corresponding to the available observations would have to be $M \times (M+n-1)$, which is always underdetermined (except when $n \equiv 1$, in which you don't have convolution of your signal).

In most convolution/deconvolution examples, the matrix is often represented as being "square," in which there exist rows in the matrix where the number of non-zero elements are less than $n$. This does not correspond to any possible measurement. Padding is common for the input signal ($\mathbf{x}$) but not for the output ($\mathbf{y}$), so how can the system not be underdetermined?

Edit for clarification:

I understand that convolution matrix is traditionally written as $N+m+1$:

\begin{equation}
  \mathbf{A} = \begin{pmatrix}
    1 & 0 & 0 & 0 & 0 & 0 & 0\\
    1 & 1 & 0 & 0 & 0 & 0 & 0\\
    1 & 1 & 1 & 0 & 0 & 0 & 0\\
    0 & 1 & 1 & 1 & 0 & 0 & 0\\
    0 & 0 & 1 & 1 & 1 & 0 & 0\\
    0 & 0 & 0 & 1 & 1 & 1 & 0\\
    0 & 0 & 0 & 0 & 1 & 1 & 1\\
    0 & 0 & 0 & 0 & 0 & 1 & 1\\
    0 & 0 & 0 & 0 & 0 & 0 & 1
  \end{pmatrix}
\end{equation}

In practice, you never measure the signals corresponding to the first two and last two rows, so it would seem that you cannot use this matrix for creating the deconvolution filter.

On the contrary, if you compute the convolution of a signal $x$ having a span of $M$ non-zero entries with a filter with a span of $N$ non-zero entries, then the resulting sequence will have $M+N-1$ non-zero entries, so that your system should be over-determined.

Computing the minimizer $x$ of $\|[0,1,0]*y-[1,1,1]*x\|_2$ leads by standard variational methods to the system of equations

$$[0,1,1,1,0]*y=[1,2,3,2,1]*x$$

with the first two and the last two removed, and for some reason mathematics / linear algebra claims that this system is solvable even if $[0,1,0]*y=[1,1,1]*x$ is not.


Adaptive Filter Gradient Descent The quadratic performance surface of an adaptive filter is a paraboloid. Its minimum can be found wherever the gradient is zero. However, since there are two types of paraboloids (elliptical and hyperbolic), is there a way to tell if the minimum detected is a global minimum or just a saddle point?

The quadratic surface is determined by the autocorrelation matrix of the data, which is always positive definite or positive semi-definite. This means that any stationary point is always a minimum. In the worst case, this minimum is not unique if the matrix is singular, but it can never be a saddle point.


Averaging over bins in frequency domain I have a sound file with a computed spectrogram from 0 to 8000 Hz (spanning about 5 octaves). The spectrogram has 128 logarithmically spaced frequency bins. How can I reduce the number of bins so that I keep an equal bandwidth in octaves? Can I average adjacent bins?

As an example, let's say I have 128 bandpass filters which are spaced equally along a log frequency axis, with center frequencies 180-7040 Hz spanning 5.3 octaves. How can I reduce this to 42 bins with equal bandwidth in octaves? Matlab or python examples are welcome.

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


Lag in financial trading vs Group Delay in signal processing I'm a newbie to Signal Processing - and so I'm providing a "background" and a "question" section separately.

Background:

Suppose I have a linear filter (w.r.t. given as

$\theta_{t+1} = \theta_t + \lambda (x_t - \theta_t)$

In Financial Trading, such a filter is often called an exponential moving average, often applied to both prices (to get smoothened prices) and to square of price returns (to get "EWMA" volatility). By conventional definition:

$\lambda =  \frac{2}{n+1}$, where $n$ is periods (example $20$-day EWMA volatility).

Usually, a lot of folks in Financial Trading define "lag" as the expected value of the filter, when $x_u$ is replace by $t-u$. For instance, the above filter can also be written as

$\theta_{t+1} = \lambda (x_t+(1-\lambda) x_{t-1}+ (1-\lambda)^2 x_{t-2}+ ...)$

To get the "lag" of this filter in Financial Trading, we write

$L(\theta) = \lambda (0+(1-\lambda) 1+ (1-\lambda)^2 2+ ...)$

which gives

$L(\theta) = (1-\lambda)/\lambda = \frac{n-1}{2}$.

Question:

Now suppose I have a filter that is non-linear in $\theta$. Example, suppose it is

$\theta_{t+1} = \theta_t + \lambda \frac{1}{2}(x_{t}^2 e^{2\theta_t}-1)$

I have no idea what would be the equivalent of a lag here. Some folks told me that the definition of lag in Financial Trading is equivalent to the group delay (or may be phase delay) in Signal Processing.

(1) I looked up the definition of group delay and phase delay (e.g. on Wiki), but I couldn't exactly relate them to the definition of lag in Financial Trading. For the linear example:

$\theta_{t+1} = \theta_t + \lambda (x_t - \theta_t)$

what would be the group delay (or may be phase delay if that's more equivalent)?

(2) For the non-linear example:

$\theta_{t+1} = \theta_t + \lambda \frac{1}{2}(x_{t}^2 e^{2\theta_t}-1)$

what would be the group delay / phase delay or the equivalent of lag as applied to the linear filters in Financial Trading?

i'm gonna work on this using the EE DSP notation i think is most conventional in the EE DSP field.  your exponentially-decaying 1-pole filter has this simple difference equation:

$$ y[n] = (1-p) \, x[n] \ + \ p \, y[n-1] $$

turns out (as we shall see below) that $p$ is the pole.  applying the Z transform gets you

$$ Y(z) = (1-p) \, X(z) \ + \ p \, z^{-1}Y(z) $$

from which we can solve for the transfer function:

$$ H(z) \triangleq \frac{Y(z)}{X(z)} = \frac{1-p}{1-p\,z^{-1}} = (1-p)\frac{z}{z-p} $$

so there is a zero at $z=0$ and a (real) pole at $z=p$.  the zero would go away if we took our input from $x[n-1]$ instead of $x[n]$ which would create a pole at $z=0$ and cancel the zero, but all that does is introduce another sample period of delay.  i don't think you wanna do that.

the (Kronecker) impulse response is

$$ \begin{align}
h[n] & = (1-p) p^n u[n] \\
& = (1-p) e^{\log(p) \, n} u[n] \\
& = (1-p) e^{n / \tau} u[n] \\
\end{align} $$

where u[n] is the Heaviside unit step function

$$ u[n] \triangleq
\begin{cases}
1, & \text{if }n \ge 0 \\
0, & \text{if }n < 0
\end{cases} $$

and $\tau$ is the "time constant" of the decaying exponential.  in order for the filter to be stable, the exponential must be decaying which means that $\tau > 0$, $\log(p) < 0$, and $0 < p < 1$

okay, so now the complex frequency response of $H(z)$ is

$$ \begin{align}
H(z)\bigg|_{z=e^{j\omega}} = H(e^{j\omega}) & = \frac{1-p}{1-p \, e^{-j\omega}} \\
& = \frac{1-p}{1-p \, \left(\cos(\omega) + j\sin(\omega) \right)} \\
\end{align} $$

the magnitude response is

$$ \begin{align}
|H(e^{j\omega})| & = \left|\frac{1-p}{1-p \, e^{-j\omega}} \right| \\
& = \left| \frac{1-p}{1-p \, \left(\cos(\omega) + j\sin(\omega) \right)} \right| \\
& = \frac{1-p}{\sqrt{ (1-p\cos(\omega))^2 + (p\sin(\omega))^2 }} \\
& = \frac{1-p}{\sqrt{ 1 + p^2 - 2 p \cos(\omega) }} \\
\end{align} $$

this is pretty familiar.  less oft observed is the phase response

$$ \begin{align}
\phi(\omega) \triangleq \arg \{ H(e^{j\omega}) \} & = \arg \left\{ \frac{1-p}{1-p \, e^{-j\omega}} \right\} \\
& = \arg \{ 1-p \} - \arg \{1-p \, e^{-j\omega} \} \\
& = 0 - \arg \{1-p (\cos(\omega) + j \sin(\omega)) \} \\
& = - \arctan \left( \frac{p \sin(\omega)}{1-p \cos(\omega)} \right) \\
\end{align} $$

the phase delay is

$$ \tau_\phi(\omega) \triangleq -\frac{\phi(\omega)}{\omega} $$

and the group delay is

$$ \tau_g(\omega) \triangleq -\frac{d\phi(\omega)}{d\omega} $$

gotta go.  i'll post this now and finish it later.  perhaps someone else wants to finish it.


with SciPy, how to generate waveform with common pulse shapes? With Python/Scipy, if I have a simply BPSK signal sequence [1, -1, -1, 1, ...], how can I generate a real-valued baseband waveforms with oversampling rate m, and common pulse shapes such as root raised cosine, Nyquist, etc.?

In Matlab I found this fdesign.pulseshaping function, is there an equivalent function in SciPy? If not how can I implement this myself?

I think in other words, how can I find filter coefficient b and a of root raised cosine or Nyquist filters so I use signal.filtfilt?

You have to write your own method to generate the root raised cosine pulse shape. This link might help you: http://en.wikipedia.org/wiki/Raised-cosine_filter

you have to select the parameters for creating one pulse shape and then you can create pulse stream for your base band signal.


Algorithm for 1d spline interpolation suitable for 8 bit microcontroler What is a concise, fast, down to earth algorithm for doing (or closely approximating) spline interpolation on a 1d continuous stream of data?
(Edit1: The paragraph below equates to saying "the data is uniform in interval." in an awkward way.)
The data is 1d in that sampling on x is at fixed, regular intervals (likely a power of 2 constant) known well ahead of time. This makes y the only variant, which should allow for quite a bit of simplification and precomputation (LUT?).
Here is a graphical example of approximately what we're trying to do. It's our algo applied to a coarsely sampled sine function.
(Edit2: Note that this is only an example of what our algo should do with coarsely sampled sine data, however, the data we would like to process will be random in the set [0,255], most likely 3 or 4 points)

(source: mathworks.com)
Assume high accuracy is not required, but that we must compute the number of results (red dots) between samples (blue circles) in less than 8ms (~120Hz). Also assume the signal processor available is limited in power and instruction set (PIC/AVR), so has only the following relevant instructions and limitations.

(signed + carry) Addition/subtraction instructions.
(unsigned 8x8-to-16) Multiplication instruction. (PIC18, megaAVR only)
Byte wide Boolean instructions (XOR, OR, NOR, AND, NAND, etc.)
Single bit left and right logical shifts. (no barrel shifter)
Can only execute at 2~4 MIPS

Additional Notes:

Responses would preferentially be in pseudocode, so they are more generally applicable.
Cheating is totally OK; it doesn't need to be perfect, just better than liner interpolation.
Bonus points for alternatives that don't need multiplication.
More bonus points for alternatives likely to complete in less than 1ms!

This is for an RGB mood lamp x-mas present for my sister and mommy :3, which I would do myself, but the maths for this are apparently beyond me.

Edit 12-21-2016: Better list formatting

Take a look at the cubic Hermite spline.  The interpolated function is continuous at the data points and the first derivative is also continuous.  Away from the data points all of the derivatives are continuous.

Let's say that the function $f(x)$ is defined by equally-spaced data points for all $x$ that is an integer.  This means you know the values of $f(0), f(1), f(2), ...$

Then separate $x$ into an integer and fractional parts:

$$ x \triangleq n+u $$

where

$$ n = \lfloor x \rfloor = \operatorname{floor}(x) $$

and

$$ u = x - n \quad \text{  ,  } \quad 0 \le u < 1  $$

$$ $$

$$\begin{align}
f(n+u) & =
\begin{bmatrix}
1 & u & u^2 & u^3 \\
\end{bmatrix}
\cdot
\begin{bmatrix}
0 & 1 & 0 & 0 \\
-\tfrac12 & 0 & \tfrac12 & 0 \\
1 & -\tfrac52 & 2 & -\tfrac12 \\
-\tfrac12 & \tfrac32 & -\tfrac32 & \tfrac12 \\
\end{bmatrix}
\cdot
\begin{bmatrix} f(n-1) \\ f(n) \\ f(n+1) \\ f(n+2) \end{bmatrix} \\
\\
& =  \frac12 \begin{bmatrix} -u^3 +2u^2 - u \\ 3u^3 - 5u^2 + 2 \\ -3u^3 + 4u^2 + u \\ u^3 - u^2 \end{bmatrix}^T \cdot \begin{bmatrix} f(n-1)\\f(n)\\f(n+1)\\f(n+2) \end{bmatrix} \\
\\
& = \frac12 \begin{bmatrix} u ((2-u) u-1) \\ u^2 (3 u-5)+2 \\ u ((4-3 u) u+1) \\ u^2 (u-1) \end{bmatrix}^T \cdot \begin{bmatrix} f(n-1)\\f(n)\\f(n+1)\\f(n+2) \end{bmatrix} \\
\\
& = \tfrac12 \bigg( (u^2(2-u)-u)f(n-1) \ + \ (u^2(3u-5)+2)f(n) \\
& \quad \quad \quad \quad + \ (u^2(4-3u)+u)f(n+1) \ + \ u^2(u-1)f(n+2) \bigg) \\
\end{align}$$

Is this enough information for how to do this in your PIC?  You need to be able to separate into integer and fractional parts and you need to be able to multiply.

in my opinion, Olli's method [now moved into its own answer] is not the best way of looking at it for the OP's case of implementing this simply in a PIC.  [his formulation] separates the four data points and computes four coefficients that are attached to powers of $u$.  that's the way to do it if your fractional ordinate is any arbitrary value that is $0 \le u < 1$.  but the OP has only a few values like $u=0, \tfrac14, \tfrac12, \tfrac34$.  or maybe 8 multiples of $\tfrac18$.

so my recommendation is to compute the vaules of these four polynomials:

$$  c_{-1} = \tfrac12 (-u^3 +2u^2 - u) \\
    c_0 = \tfrac12 (3u^3 - 5u^2 + 2) \\
    c_1 = \tfrac12 (-3u^3 + 4u^2 + u) \\
    c_2 = \tfrac12 (u^3 - u^2)  $$

and do that for every fractional value of $u$ (such as $u=0, \tfrac14, \tfrac12, \tfrac34$) that you will use many, many times.

then the code in the PIC only needs to implement is a dot product between the 4 data points and the selected set of coefficients:

$$ f(x) = f(n+u) = c_{-1} f(n-1) + c_0 f(n) + c_1 f(n+1) + c_2 f(n+2) $$

since $c_{-1}$ and $c_2$ can be shown to always be negative for $0 < u < 1$, then put into the table their absolute values and subtract their terms:

$$ f(x) = f(n+u) = c_0 f(n) + c_1 f(n+1) - (-c_2) f(n+2) - (-c_{-1}) f(n-1) $$

the coefficients stored will be 256 times bigger than their actual value (and stored as 8-bit unsigned integers) then, after multiplying and accumulating your answer (that is 256 times too big), you add 128 (for rounding) and shift right 8 bits (which is the same as taking the answer out of the higher-order byte of the product).


Can a narrowband beamformer be used for speech separation In my master thesis I am trying to implement a beamformer.

I've already done the delay-sum part 1 year ago. I have been told:

in FIR filtering we can do

FFT --> null the unwanted frequency bin --> IFFT .

We can to the same in beamformers. If we know the interference DOA (like 0 in angle 25 and 1 in everywhere else between -90 to 90) output we can ---> IFFT and use the weights in aparture elements.

Whatever I tried it didn't work. Than I begin to read books , surprisingly couldn't find the method that my advisor told me. Instead I found out methods like MVDR, LCMV. Now I am tring to use matlab method lcmvweights to get correct weights for each element and apply in my delay-sum beamformer.   Even I succeed to use this method I am curious if narrowband-beamformers could be used in complex signals like speech?

If you are dealing with say 8 kHz for a nominal speech bandwidth of of 4kHz i.e. 0 Hz - 4 kHz, then the speech is essentially a wideband signal. Therefore narrowband beamforming won't work very well. Your beam pattern will be okay for the particular frequency of interest - but once you start moving away from that frequency your beam patterns will deteriorate.

What you are looking for is wideband beamforming - which uses actual time delays (or linear phase shifts across frequency in the frequency domain) rather than just phase shifts.

Usually the techniques like LCMV and MVDR are developed for narrowband signals. There are a couple of ways of approaching the wideband problem:

Use a tapped delay line on each channel. If you have $n$ channels and $m$ taps per channel then your correlation matrix is $mn$ x $mn$. So the systems of equations gets very large.
Use a bunch of narrowband beamformers. In this case you would have $m$ beamformers ($m$ frequencies) each one having $n$ channels. Now each beamformer has a $n$ x $n$ correlation matrix, but you have $m$ of them. So it leads to a reduction in complexity from the previous case.
You can form a set of $b$ conventional beams (using time-delay rather than just phase multiplier) and then do the adaptive processing on the beams.

The best reference I can think of is by Van Trees - Optimum Array Processing. Note you sometimes run into slight differences in terminology - some texts will denote MVDR as using the signal correlation matrix, while others use the signal plus noise correlation matrix - Just be careful which one you're looking at. I know Van Trees makes a distinction between the two cases. Most other texts do not and only use one formulation and call it MVDR.


Can $\delta(t+\infty)$ be a legitimate signal? Mathematically speaking, when I try to use some signal to disprove a system is invertible, can I use the signal like $\delta(t+\infty)$ ($\delta$ representing the Dirac distribution)? For example, the two inputs $0$ and $\delta(t+\infty)$ are distinct, but their outputs from the system $\int_{-\infty}^{t}e^{\tau}x(\tau)d\tau$ will then both be $0$. In this way, can I say the system is invertible?

It's possible, in mathematics, to complete real numbers with "infinite" values, with sound topological properties; for instance non-standard analysis or the Extended real number line (discussion at Math StackExchange).

However, in this context, for any standard real number $t$, the rule is to set $t\pm \infty = \pm \infty$ (see Arithmetic operations). So, using your notations, for every $t$, you would have [CAUTION ADVISED] $\delta(t+\infty)=\delta(\infty)$, to which I am not able to give a meaning, else than something being a constant evererywhere equal to $\delta(\infty)$. Of course, $\delta$ is a distribution, and should not be treated as a function, but defined as a functional operator, or via limits of functions.

However, so far, I have never seen

a (serious and useful) use of the extended real system in signal processing,
a reference trying to define such a Dirac at infinity.

For the latter, I fear (just intuition) that constructing this with test functions can be troublesome, since two different limits are required (location and amplitude). This is just an analogy, but $\delta(0)\times \delta(0)$ is, as far as I know, not defined.


Linear Combination of Two Signals In my signal processing class I have some basic tasks working with sinus signals in Matlab. I have two sinusoidal signals with 1kHz and 3kHz and sampling rate of 16kHz. The next task says 'Create a linear combination of the two signals and save the result in a new vector'.

I don't really understand what's meant by the linear combination. I remember something vague about having to add them with weights but I'm don't quite comprehend what that does or how to do it.

As you said, a linear combination of two signals $x_1[n]$ and $x_2[n]$ is simply

$$y[n]=w_1x_1[n]+w_2x_2[n]$$

with real-valued weights $w_1$ and $w_2$.


Adding echo to audio sample I have recorded my own voice in Matlab and I intend to add some echo to it.I came up with one solution for getting the desired echo effect:

Delay the sampled audio in the time domain and adding it to the original sample.

In order to do this I'm using Matlab and I have basically done the following:

recObj = audiorecorder(44100, 16, 2);%sampling rate of 44100Hz, stereo
recordblocking(recObj,length);%record audio for a fixed length duration
y = time_delay(getaudiodata(recObj) , 5000 );%set a delay of 5000 to original sampled audio data

function [ y ]  = time_delay ( x , R )
%this function sets a delay to x of R
y = zeros(length(x) + R , 1);
y(R + 1:length(y)) = x(1:length(x));

As you can see I have a recObj wich I can use for playing the recorded audio.In order to get a delayed version of the recorded audio I have used getaudiodata(recObj) and my own function time_delay.A problem that I have is that I don't know how to get a new delayed audiorecorder object from the delayed audio data retrieved by the time_delay function and I need an audiorecord object cause this is the only way of playing something.
So, my question is:

Is this a good approach to get the desired echo effect in an audio
signal?
If it is, then how do I set the delayed audio data to an
audiorecorder object?

Here is the code for adding the two signals(the delayed and not-delayed):

    x = getaudiodata(recObj); n1 = 1:size(x,1);%audiodata of original signal
    y = time_delay(x , 50000 ); n2 = 1:size(y,1);%audiodata of delayed signal
    mixed = sigadd(x,y,n1,n2); %audiodata of mixed signal
    mixrecObj = audioplayer(mixed,44100);
    play(mixrecObj);

This is the sigadd function:

function [ y,n ] = sigadd( x1,x2,n1,n2 )
%implements y[n] = x1[n] + x2[n]
%y = sum sequence over n wich includes n1 and n2
%x1 = first sequence over n1
%x2 = second sequence over n2
%

n = min(min(n1),min(n2)):max(max(n1),max(n2)); %duration of y(n)
y1 = zeros(1,length(n)); y2 = y1;
y1(find((n >= min(n1))&(n <= max(n1)) == 1)) = x1;
y2(find((n >= min(n2))&(n <= max(n2)) == 1)) = x2;
y = y1 + y2;
end

This works perfectly, but is not enough for getting a real echo effect. In order to archieve this we have to
create more delayed signals(with diferent delays of course), and then add them all,like this:

x = getaudiodata(recObj); n1 = 1:size(x,1);
y = tim_delay(x , 5000 ); n2 = 1:size(y,1);
s = time_delay(x , 4000 ); n3 = 1:size(s,1);
d = time_delay(x , 3000 ); n4 = 1:size(d,1);
mixed1 = sigadd(s,d,n3,n4);n5 = 1:size(mixed1,2);
mixed2 = sigadd(x,y,n1,n2);n6 = 1:size(mixed2,2);
totmixed = sigadd(mixed1,mixed2,n5,n6);
mixrecObj = audioplayer(totmixed,44100);
play(mixrecObj);

I have tested this one with the three delays plus the original and you can hear a pretty nice echo.

Note: In the third piece of code it would have be better if there were a function that returned the signal with the echo by passing an echo effect index  as a parameter to it, but I left it that way for the sake of simplicity.


Fourier transform for non-sinusoidal signal Using MATLAB's fft function, I would like to retrieve the amplitude and frequency of this particular signal :

My data file (a) has 2 columns, first one for time ans the second for amplitude, I tried this script

    DT = (a(2,1)-a(1,1));
    Fs = 1/DT; % sampling frequency
    NFFT = 2^nextpow2(n);
    xdft = fft(a(:,2),NFFT)/n;
    f = Fs/2*linspace(0,1,NFFT/2+1);

    plot(f,2*abs(xdft(1:NFFT/2+1)))

The output is :

which is not corresponding to my amplitude. Thanks a lot for your help.

fft with all the harmonics :

The Fourier Transform will decompose your non-sinusoidal signal into harmonics, dominantly odd harmonics since your distortion appears symmetrical, and the amplitude as you derive would be the amplitude of the relevant sinusoidal harmonic (so in your case it looks like the fundamental is shown, so we are seeing the amplitude of the first fundamental harmonic which is a sinusoidal signal). It will be less than the peak amplitude shown in the time domain plot, which is the composite of all the harmonics.

Other factors that will effect the amplitude is spectral leakage due to your use of a rectangular window (so to the extent the harmonics fall into a sidelobe of the kernel for your rectangular window) and scalloping loss (to the extent the fundamental frequency is between an integer sub-multiple of your sampling frequency. Both of these effects are described in more detail in my favorite paper by fred harris: fred harris On the Use of Windowing

I am surprised that we do not see any evidence of these harmonics in your frequency plot but this may be because you are only showing a portion of all the frequencies or perhaps because the magnitude is not on a log scale we are not able to make out the harmonics. For example, there should definitely be a signal at three times the fundamental shown. The fundamental appears to be close to 0.04; meaning we should see a harmonic at 0.12; if your frequency of 0.1 represents the Nyquist boundary, then the image of this would be at 0.08 assuming a real signal as you have plotted. (Could you please udpdate your plot to have a log (dB) magnitude scale so we can be sure nothing else is astray?).

That said, we do see that your amplitudes are actually quite close visually from the plot (approx 0.17 in the FFT vs approx 0.19 in the time domain plot), so you must be concerned with the small difference, which is accounted for from the effects described above (the 0.17 is the amplitude of the primary tone as modified by the spectral leakage of your rectangular window and any scalloping loss).


Application of notch filter in radio receivers So according what I have read so far notch filter (or stop-band filter) can be used to get rid of the narrowband noise when it is present in the frequency spectrum of the received signal.

For example stray sine waves may impinge on the receiver as we listen to radio.

My question is, how would the notch filter detect the particular noise frequency in first place so that it could remove it?

For a given system and application, especially wideband sytems, sometimes we can consider dominant interfering signals. So a notch filter is tuned specefically on the frequency of the interferers. In such case, we remove part of the desired signal's energy as well, however, the narrowband interference wich is usually much stronger than the desired signal (in a narrow bandwidth only) is removed and its harmful effect is prevented.

As an example, you might have heard of ultra widebsnd systems that span over a large spectrum between 3.1 Ghz to 10.6 Ghz. They may interfere with WLAN at 5.6 Ghz. Therefore, the UWB receiver can deploy a notch filter at 5.6 GHz.


Evaluating loudspeakers If you don't have a lot of money to spend but you want to estimate the frequency response of different loudspeaker setups what is the best way to do it?

I bought a very nice microphone and built a box with foam on the inside. My idea is to place the loudspeaker and the mic in the box and measure the power level at different frequencies by playing a sine tone at those frequencies.

Or is it better to have the microphone and loudspeaker in an open space where the mic is 10 inches away from the loudspeaker?

I'm sure neither of the 2 methods above are valid for making such measurements (if so please explain why), but how should I do then?

Also, should I use sine tones or is it better to use another test signal like an "impulse", "white" noise or a chirp signal?

Hoping to get some good answers / advice

Thank you!

Not an easy job. A loudspeaker produces a complicated 3 dimensional sound field, which heavily interacts with the environment. Even at a distance of a about 3 feet the reflected energy is significantly bigger than the direct energy at most frequencies. Typical loudspeaker measurements are

On axis response in an anechoic chamber
Full 3D response in an anechoic chamber
Radiated power in the far field in a highly reflective environment.
Room responses taken in "real world" setups, deriving a perceptually relevant spectral metric and averaging over a large number of environments and placements

Each of these measurements require specific measurement methods and techniques which are also not trivial.

Anything you can do with a single microphone in real room will most likely be either useless and misleading. The box is also not a good idea since any box has very strong internal acoustic modes and what you will be measuring will be dominated by these modes regardless of what the loudspeaker is doing.

Depending on the physical size of the loudspeaker and the size of your measurement space, you may be able to pull off a near field measurement with careful time gating of the impulse response (loosing a lot of low frequency resolution in the process). Getting this right is a lot of work and in my opinion not worth the bother.