How do I numerically calculate a function from its noisy gradient? I have the model $\ s(x,y)=x^2+y^2, 0 \leq x \leq 1, 0 \leq y \leq 1 $.

Instead of observing the model directly I am observing the derivatives of the model + some noise (e):

$\ p(x,y)=s_x+e, q(x,y)=s_y+e$

From measurements of p(x,y and q(x,y) I want to estimate s(x). Say I know that s(0,0)=0.

According to the gradient theorem:
$\ s(x,y) = \int_{(0,0)}^{(x,y)} [s_x,s_y]d\vec{r}$

regardless along which path we integrate.

As a small experiment (in Matlab) I added normal distributed noise, N(0,1), to p=2x and q=2y. Then I integrated first along x followed by along y: SXY. Next I integrated first along y followed by along x: SYX.

The results show that the gradient theorem do not hold in this case (because of the noise):

The Root Mean Square Errors relative to the model are:

ErmsXY =
    0.1125
ErmsYX =
    0.0920

How can I find a better (less RMS error and smoother) estimate of s from p and q?

EDIT:

From what I read; using the curve integral is referred to as local integration.
There are also global integration methods where one instead tries to choose a S(x,y) that minimizes:

$\ \int_0^1\int_0^1 [|S_x - P|^2  + |S_y - Q|^2] dxdy$

Global integration methods are supposed to give better results when the gradient is noisy,
but how do I do this in practice?

EDIT 2:

One approach that I have used is this:

first we introduce linear derivation operators: $\ s_x = D_x*s,s_y = D_y*s$.

The result is the following linear equation system:

$\ D_x*s=p, D_y*s=q$

Next find a Least Square Error solution to these equations.
A LSE solution to these equations is supposed to be equivalent to minimizing the integral from above.
How can this be shown?

The results are good:

The RMS error is about 1/5 of that of SXY and SYX and the solution is also smoother.

However there are some drawbacks to this approach:

it is difficult to implement; must use central differences and "flatten" 2D s matrix into vector etc.
The derivation matrices are very large and sparse so they may consume a lot of RAM.

Another approach that seems potentially both simpler to code, less RAM consuming and faster is to use FFT. In the Fourier space these pdes become an algebraic equation.
This is known as the Frankot-Chellappa algorithm, but unfortunately I have not got it to work on my example data.

You can filter either the gradient itself or the result, $s$.  You would need to know the characteristics of the true gradients fairly well to know what the frequency bandwidth is.  At that point you could design a low-pass filter that would preserve the signal but get rid of higher frequency noise.


phase tracking loop for dbpsk I have a communication with a DBPSK modulation. The communication is disrupted by a Doppler shift (due to movement) and Doppler spread (due to channel). At the reception, after Doppler shift compensation with a preamble and DBPSK demodulation, I have a problem with the phase of the symbols. To correct the phase error, I use a phase tracking loop. Here is the link on what I'm doing phase tracking loop.

Here are my questions:

1) How to choose the constant loop $K_c$?

2) We suppose the Doppler shift is perfectly compensated. In my Matlab code, if for example I take $K_c=0.01$, sometimes I find exactly the symbol sent and sometimes I have symbol errors. What is the reason for these errors?

Thank you

Ok, to get this straight: you're doing DBPSK, but you're still seeing a constant phase error that has a non-zero mean?

That means you've got a frequency error.

Your assumption

  We suppose the Doppler shift is perfectly compensated

hence must be wrong.

So, correct frequency instead of the phase. You're already using DBPSK because you don't want to correct the absolute phase (I hope this was known to you. If you didn't know why you did differential modulation: please go and read up on differential phase shift keying, and why you don't use a phase-correcting loop with D(B)PSK).


Extract signal from big gaussian noise I have signal, periodic with amplitude 10 with frequency about 10 kHz. It is hidden in Gaussian noise with standard deviation 100. Some idea how to extract periodic signal?
Thanks in advance.

Rearrange the signal into a matrix containing one period in each row and average down the columns.


Properties of Invertible System I have come to the understanding that if these three properties hold true then the system is a memory system:

1) Time Scaling:

$$y(t)=x(2t)\,,$$

2) Time Shifting:

$$y(t)=x(t+2)\,,$$

3) Integration based system.

I want to ask that do these properties also hold true for invertible systems?

In the context of signal processing, a continuous-time system

$$y(t) = \mathcal{T} \{ x(t) \}$$

is said to have memory if its output at time $t=t_0$ has a dependency on the input values for any time other than $t=t_0$. Otherwise; if the output at time $t = t_0$ depends only on the input value at time $t=t_0$, then it's a memoryless system.

So the system defined by:
$$y(t) = x(2t)$$
has memory, which can be proven by a single example, as $y(t)$ at $t=5$ is given by $x(t)$ at $t=10$...

Furthermore, a system is invertible if unique (distinct) inputs produce unique (distinct) outputs, otherwise it not invertible.

So the continuous time system defined by:
$$y(t) = x(2t)$$
is invertible, as its inverse is given by $x(t) = y(t/2)$.


System Stability: Can we derive stability of a discrete system (Frequency domain, Z-transform) by applying analogous methods? So given some analogue system function in the complex s-domain. Can we perform a stability analysis in the $s$-domain, before actually transfer it into the $z$-domain? So in other words analysis in the $z$-domain is not necessary anymore?

The reason why I am unsure is that we have a sample and hold element which needs to be added to our system:

$\frac{1-e^{-sT}}{s}$ which changes to $\frac{z-1}{z}$

(I did use partial fraction decomposition/method of residues and did not use any approx. like bilinear transform, impulse variance method.)

If you use the bilinear transform to map from s-plane poles and zeros to z-plane poles and zeros, the left-hand half-plane of the s-plane will be mapped to the inside of the unit circle on the z-plane. This preserves stability because on s-plane the requirement for stability is that poles are on the left-hand half-plane and on the z-plane it is that they are inside the unit circle.

I did not examine too carefully your s and z-domain formulas, but indefinite integration is marginally stable with a pole at the origin of the s-plane, or on the unit circle of the z-plane. Implemented as a nonrecursive filter:

$$\text{out}[i] = \text{in}[i-N + 1] + \text{in}[i-N + 2] + \dots + \text{in}[i-1] + \text{in}[i],$$

integration over $N$ samples is stable. Implemented as a recursive filter:

$$\text{out}[i] = \text{out}[i-1] + \text{in}[i] - \text{in}[i-N],$$

it is stable in fixed-point arithmetic, but can drift if implemented in floating point arithmetic due to accumulation of numerical error. Even then it can be stabilized by moving the pole slightly inside the unit circle, or by periodically resetting the filter state by calculating it using the nonrecursive filter. But sample and hold is about taking an input sample and keeping the output at that value for some time. That can be programmed directly and gives no stability problems.


Image Denoising with Better Edge Preservation I have the input image :

and the output of vein detection for the leaf using a Gabor filter, but the output is really noisy:

I tried using Total variation denoising however the results are not good:

However I don't want to loose the fine details in the leaf's veins, so a median filter won't suit my problem

Sounds like you want to denoise and preserve edges. Have you considered nonlocal means? There's some GPL'd C++ code along with a brief writeup of the algorithm by the original authors here: http://www.ipol.im/pub/algo/bcm_non_local_means_denoising/

One caveat, nonlocal means is very slow and the output can be sensitive to the implementation you have. You may also consider ROF minimzation as it's fast and does a good job of preserving edges. Here's some matlab code that does it: http://www.stanford.edu/~tagoldst/Tom_Goldstein/Split_Bregman.html


Signals cross-correlation I am trying to measure phase difference between two received waves which are recorded by two separate PC's (microphones). The two PC's starts recording at the same time, then I make a sine wave sound with frequency of $1000\textrm{ Hz}$ using a sound speaker which is recorded by the two microphones.

For the purpose of measuring the phase difference between the two recordings, I do a cross-correlation between the two waves. Theoretically ... the distance between zero and the index of the peak value of the output graph (after cross-correlation) is the phase difference between the two waves.

For some reason I am getting a graph with more than one peak. I think this is because of the noise.

Does someone know if I can make a normalized cross-correlation to sound waves, does this solve my problem  ?
If yes, how can I do this ?
If no, any suggestions to solve this problem ?

The waves and the correlation after plotting in MATLAB:

I know the two waves are not perfectly the same, they represent a sine wave with $1000\textrm{ Hz}$, which is recorded by two separate desktop PC's.

As you can see, there are 3 peaks in the cross-correlation graph. And this is just an example , some recordings gave me about 10 different peaks.

When using a constant tone audio beacon, beware of room echoes causing multi-path interference and distortion, especially around the leading and trailing portions of your received waveforms.

Try using a frequency sweep instead of a constant tone for your transmit waveform.  This might provide you with a sharper correlation peak that is less likely to have duplicate peaks.


Signal normalization I have some questions about what is normalization. The obvious way to do that is

x_max = max([max(x), abs(min(x))]);
x=x-mean(x);
x_norm=x./x_max;

Can somebody explain how it is done?
And could you please explain what is the benefit of this method?

[EDIT] After a second read, the proposed normalization looks non standard. Suppose that $m\le x \le M$ ($m$ and $M$ denote the min and max). The scaling factor will be, depending on the situation:

if $m\le M\le 0$: $-m$,
if $m\le 0\le M$, and $|m|\ge M$: $-m$
if $m\le 0\le M$, and $|m|\le M$: $M$
if $0\le m\le M$: $M$

It turns out to be (if I do not err) the largest of the absolute values of $m$ or $M$.
A more standard writing, more symmetric, could be:

x_max = max([abs(max(x)), abs(min(x))]);

or directly:

x_max = max(abs(x));

This computes the maximum x_max in absolute amplitude. This corresponds to an $\ell_\infty$ norm dispersion-type. Then you center $x$ around the mean (in some $\ell_2$ sense), then divide by the above maximum x_max.

If $\mu$ denotes the mean, you will end up:

for the first two cases, in $[\mu/m-1,\mu/m-M/m]$
for the last two cases, in $[m/M-\mu/M,1-\mu/M]$

And there is not a lot more I can say, without more information. If the extrema values are sufficiently symmetric around the mean, with a close-to-$0$ average, the outcomes would lay  around the $[-1,1]$ interval, because $M/m \approx -1 $ and $\mu/M$ or $\mu/m$ are close to $0$. But you could have different behaviors: if $m=-10$, $M=1$, $\mu=-9$, you will have x_norm in $[-0.1,1]$.
But if the mean is $0.5$, the interval would be $[-1.05,0.05]$.

Variable reduction  is something one could take classically as an input to a decorrelating or source separation technique;  a PCA for instance, or related tools (SVD, independent component analysis) for instance.
But subtracting the central location (here the average) before computing extremas and rescaling a dispersion after that is more standard.

And however, mixing different norms in normalization is not the practice I would recommend. So:

subtracting the mean and dividing by the (residual) standard deviation
subtracting the median and dividing by the MAD (median absolute deviation)
subtracting the midrange and dividing by the max or min of the residual (ie the half range)
mapping the values between fixed bounds, such as $[0,1]$

are much more common.


How can create a .wav bar based spectrogram on matlab? How can I create a bar based spectrogram that updates as a .wav file is playing. This is commonly found in music players. Here is an example below:

I am able to display a amplitude/frequency graph that updates as the .wav file is being played using this code below:

function syncPlayerDemo()
%some example music
%set up audio player
%player = audioread('test.wav');
[y,Fs] = audioread('test.wav');
player = audioplayer(y, Fs);
[samples,channels]=size(y);
%calculate timeline
t=linspace(0,1/Fs*(samples-1),samples);
%initialize full plot, update will only move the visible area using xlim
h=plot(t,y);
%set up callback to update every <TimerPeriod> s
player.TimerFcn=@timerFcn;
player.TimerPeriod=0.1;
player.playblocking()
end
function timerFcn(source,data)
%an area of length <area> s will be visible
area=1;
position=(source.CurrentSample-1)/source.SampleRate;
%move visible area, current position is in the center
set(gca,'XLim',[position-area/2,position+area/2]);
%used a waitbar for testing, might be commented in
%waitbar(source.CurrentSample/source.TotalSamples);
end

Which results in an updating graph:

How can I get the bars like above, any suggestions?

You can get the bar based spectrogram by decomposing the signal segments into fractional-octave subbands. Try octaveFilterBank and then compute signal power level by pow2db(sum(abs(signal))). This will result in a logarithmic frequency scale.
You can also use pwelch to estimate power spectral density, which brings you a linear frequency scale.


Generalized Cross Correlation What is the definition of Generalized Cross Correlation?

I need it for estimating the delay (time offset) between two audio signals. The time offset between the two signals might be as large as 500ms.

EDIT  The OP mentions below Deve's answer that this was not the "Generalized Cross Correlation" algorithm referred to.  The real one seems to be this one:

This paper seems to give equations about how to implement it in the time domain.

The equation they give is basically what Deve reports below. The key point they make is that the "input" to the cross-correlation is not the bare signals, but the residuals of an LMS algorithm applied to the bare signals.

This as the effect of pre-whitening the input signals.

The method seems to be from this paper:

Azaria, M.; Hertz, D., "Time delay estimation by generalized cross correlation methods," Acoustics, Speech and Signal Processing, IEEE Transactions on , vol.32, no.2, pp.280,285, Apr 1984

doi: 10.1109/TASSP.1984.1164314

Abstract: The problem of estimating time delay by cross correlation methods is reexamined for the whole class of stationary signals. Expressions are derived for the estimation mean square error (MSE) by the cross correlation method, and are shown to be identical to previously published results for Gaussian signals. The generalized cross correlation method is also analyzed, and the optimal weight function for this method is derived. It is shown to be identical to that derived for Gaussian signals by the maximum likelihood method. For the cross correlation method a simplified MSE expression is derived, which is to be used instead of a previously published result.

keywords: {Autocorrelation;Bandwidth;Correlation;Delay effects;Delay estimation;Error analysis;Estimation error;Maximum likelihood estimation;Mean square error methods;Signal to noise ratio}

From which the key formula appears to be:

where


Preprocess sound for voice recognition in Pocket Sphinx on Raspberry Pi 3 I want to improve my voice recognition system on a Raspberry Pi 3 running Pocket Spinx. It's pretty accurate (limited recognition vocabulary) but is terrible when other noise is around.

Is there any way to pre-process sounds before feeding them into pocket sphinx to filter only the human voice ? I was thinking there might be some way to reduce the data being passed in, eg human voice being between certain frequency etc, white noise removal etc. I haven't found anything simple so far. I'm thinking of filtering and processing the voice recordings on an STM32 that has an Arm Cortex M4 and DSP etc.

Human voice has the same frequencies as other noises, it is not possible to filter it with dumb algorithm. It is possible to filter based on temporal patterns like spectral subtraction, but this way you can only filter static noise, not background speech/bable noise. Pocketsphinx already has spectral subtraction.

To filter noise reliably you need an advanced noise filtering neural network or something like that, see RNNNoise, it is very resource-intensive. Also, it is better to recognize noise sound instead of filtering it since recognizer can figure out things by himself.

One important approach is to reliably filter is to use direction of arrival of the sound, that is a microphone array which can filter noises based on direction. An example of such array is Respeaker or Matrix One. Thats why all smart speakers have microphone arrays as a critical component.

For more accurate recognition it is better to use modern toolkit with neural network like Kaldi, it is also much more accurate with noise filtering and without it.

Further reading for you:

Snips - Snips is a partially proprietary framework, but you can get some nice ideas and hardware advices. It is easy to get a good system using their packages too. Also see Community projects
Make a smart speaker with Respeaker
Zamia Speech the easiest way to run Kaldi on Raspberry Pi


Cascading a filter with itself to form a higher-order filter? I've been given an audio signal with some high frequency noise. The professor gave me the system function of a lowpass filter that gives me control over the cutoff frequency:

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

where $\alpha=\frac{1-\sin \omega_c}{\cos \omega_c}$ and $\omega_c$ is the cutoff frequency in radians/sec.

The professor claims that if the filter is cascaded with itself $L$ times, it becomes a filter of order $L$.

My questions are the following:

What does increasing the order of the filter do?
Why would I want to increase the order of the filter?
If I cascade the filter with itself $L$ times, won't I lose the ability to control the cutoff frequency? How could I correct for this?

edit:

Here's what I get for the frequency response of the filter when the cutoff frequency is set to 1000 Hz. Sampling frequency is 16000 Hz.

Replace your $\alpha$ with the following to see if it helps:

fs = 16000;
fc = 3500;
wc = 2*pi*fc/fs;

L = 5;            % order of the cascade
B = 2^((L-1)/L);
C = cos(wc);
M = B-C-1;
N = 1+C-B*C;
alfa = ( -N + sqrt( N*N - M*M ) )/M;

I've got the following plots, compared to your posted $\alpha$ calculation.


How to merge transitions? I need to represent my signal as a set of transitions from one level to another, and I am interested in transitions that have relatively big amplitude 2000+.

I applied median filter to get rid of the high frequency part of the signal and result looks like this. Now I can take gradient and get set of transitions, and based on the magnitude of these transitions, keep only big one. So, I will get transition, corresponding to the left side of this peak, I will filter out transitions from the top. But I would like to have right side of this peak also to be represented as one nice 2000+ transition. But as you can see on the picture, there is small kink in the middle of the right side that will break my transition into 2 parts.

Is there some kind of a filter that will transform my signal and remove this kink?

One of many possible algorithms to provide such edge sensitive information to its users are classified under change-detectors. The one that best fits to your kind of data uses thresholds, accumulators, memory, leakage and smoothing factors as well as variance and mean normalizers to decide on which changes belong to meaningful edges and which are just local variations.


Why do we choose 44.1 kHz as recording sampling rate? Peoples' ears can hear sound whose frequencies range from 20 Hz to 20 kHz. Based on the Nyquist theorem, the recording rate should be at least 40 kHz. Is it the reason for choosing 44.1 kHz?

It is true that, like any convention, the choice of 44.1 kHz is sort of a historical accident. There are a few other historical reasons.

Of course, the sampling rate must exceed 40 kHz if you want high quality audio with a bandwidth of 20 kHz.

There was discussion of making it 48.0 kHz (it was nicely congruent with 24 frame/second films and the ostensible 30 frames/second in North American TV), but given the physical size of 120 mm, there was a limit to how much data the CD could hold, and given that an error detection and correction scheme was needed and that requires some redundancy in data, the amount of logical data the CD could store (about 700 MB) is about half of the amount of physical data. Given all of that, at the rate of 48 kHz, we were told that it could not hold all of Beethoven's 9th, but that it could hold the entire 9th on one disc at a slightly slower rate.  So 48 kHz is out.

Still, why 44.1 and not 44.0 or 45.0 kHz or some nice round number?

Then at the time, there existed a product in the late 1970s called the Sony F1 that was designed to record digital audio onto readily-available video tape (Betamax, not VHS). That was at 44.1 kHz (or more precisely 44.056 kHz).  So this would make it easy to transfer recordings, without resampling and interpolation, from the F1 to CD or in the other direction.

My understanding of how it gets there is that the horizontal scan rate of NTSC TV was 15.750 kHz and 44.1 kHz is exactly 2.8 times that. I'm not entirely sure, but I believe what that means is that you can have three stereo sample pairs per horizontal line, and for every 5 lines, where you would normally have 15 samples, there are 14 samples plus one additional sample for some parity check or redundancy in the F1.  14 samples for 5 lines is the same as 2.8 samples per horizontal line and with 15,750 lines per second, that comes out to be 44,100 samples per second.

Now, since color TV was introduced, they had to bump down slightly the horizontal line rate to 15734 lines per second. That adjustment leads to the 44,056 samples per second in the Sony F1.


Oversampling then filtering and downsampling Currently I'm working on a project to record ELF (Extremely Low Frequency) signals from a receiver. The highest frequency of interest is about 30 Hz and the system should have a flat frequency response.  I'm looking for a simple solution using an 8 bit ADC and I wonder how the resolution can be improved using oversampling and to alleviate the requirements for the analog low pass filter before the ADC. Without oversampling and a sampling rate of 64 sps a rather steep analog low pass filter is needed to avoid aliasing. With oversampling and a gain of 3 more bits the sampling rate is 2^(2*3) times 64 sps = 4096 sps and the analog filter before the ADC is much easier to realize. The problem is to convert the 4096 sps signal back to 64 sps using a digital low pass filter and downsampling.

From what I read so far, the attenuation needs to be almost 68 dB at 32 Hz for 11 bits but at the same time signals up to about 30 Hz should pass. It seems a very steep digital low pass filter is needed. I tried to do some calculations with filter design software, but the required order is too high to get results. Is this really so hard to realize or am I missing something? Any help would be greatly appreciated.

As you've figured out, for the filter specifications you gave, you need a very high filter. To summarize, you're asking for:

Passband edge: 30 Hz
Stopband edge: 32 Hz
Minimum stopband attenuation: 68 dB

The main determining factor for the filter order that you need to meet the specifications is the width of the transition band as a fraction of the sample rate. In your case, you want a transition band that is 2 Hz wide, at a sample rate of 4096 Hz. That's a very sharp transition!
A better approach is to perform the decimation in multiple filter stages that are cascaded in series. This allows you to meet the overall frequency response requirements that you have without using a single very-high-order filter. You need an overall decimation factor of 64 to be implemented by the cascade of filters. It can be hard to determine the optimal (in terms of total overall computations to apply the filters) arrangement of stages, but I would recommend the following approach as a start:

Stage 1: Design a filter with a passband edge at 30 Hz and a stopband edge of (512 - 32) = 480 Hz. Decimate the output of the filter by 8. Thus, at the output of this filter, you will have a signal sampled at 512 Hz.
You might notice that this approach yields some aliasing, as the transition band of the stage-1 filter is allowed to extend beyond the new Nyquist frequency of 256 Hz. However, this isn't a problem, because none of the aliasing occurs within the desired passband of the overall cascade. You're allowing some energy to alias through into what will become the stopband of stage 2, where it will be eliminated. Stretching out the transition band of the stage 1 filter allows you to save some computations, because the wider transition will require a lower filter order. Since stage 1 operates at the highest sample rate, this can be a significant benefit.

Stage 2: Design a filter with a passband edge at 30 Hz and a stopband edge of 32 Hz. Decimate the output of the filter by 8. Thus, at the output of this filter, you will have a signal sampled at 64 Hz with a 30 Hz passband.
While this is still quite a sharp transition band, it should be easier to design because the transition band width is a factor of 8 larger when expressed as a ratio of the filter's sample rate. In addition, the work required to apply the stage-2 filter is reduced further because it operates at the lower 512 Hz rate instead of 4096 Hz.

Using MATLAB's fdatool, the above specifications yield the following filters:

Stage 1: 29th-order FIR, 0.1 dB passband ripple, 80 dB stopband attenuation

Stage 2: 852nd-order FIR, 0.1 dB passband ripple, 80 dB stopband attenuation

I just used the default attenuation settings for illustration. You'll notice that with this setup, the second-stage filter is still quite complex. While I've used filters of this length before, it's likely that you might want to break the problem down even further to avoid such long filters. To do so, you would add more stages of decimation. For instance, you might try three stages of decimation by 4, or 6 stages of decimation by 2. The best way forward is dictated by your system's constraints, but you should see the idea of the approach.


Is an FIR filters 'transition width' frequency dependent? I was recently brushing up on some FIR filter design and was browsing the Wikipedia page. A claim in the article caught my attention:

  The main disadvantage of FIR filters is that considerably more
  computation power in a general purpose processor is required compared
  to an IIR filter with similar sharpness or selectivity, especially
  when low frequency (relative to the sample rate) cutoffs are needed.
  However many digital signal processors provide specialized hardware
  features to make FIR filters approximately as efficient as IIR for
  many applications.

Article

Now what caught my attention was the statement made about how more computation is required when a low cutoff frequency relative to the sample rate (i.e. low normalized cutoff frequency) is specified.

My question: is this true? Or is it just a poorly written statement that assumes you need sharp transitions when working with low cutoff frequencies?

I was always been under the impression that FIR filters become very difficult to implement when trying to achieve sharp transitions, not certain frequencies. Have I been wrong in this assumption?

Furthermore, this article seems to confirm the fact that the decay rate of the impulse response is independent of a chosen cutoff frequency. So assuming the impulse response has a constant decay rate, shouldn't an FIR filter of fixed length being equally good (or bad) regardless of cutoff frequency?

You're right: the required filter order is approximately inversely proportional to the desired transition bandwidth $\Delta\omega$, regardless of the cut-off frequency. This is reflected in the empirical formulas for estimating the required filter orders for the Kaiser window design method as well as for the Parks McClellan equiripple design of low pass filters:

$$\begin{align}M&=\frac{-20\log_{10}\delta-8}{2.285\Delta\omega}&\qquad\text{(Kaiser Window)}\tag{7.63}\\{}\\
M&=\frac{-10\log_{10}(\delta_1\delta_2)-13}{2.324\Delta\omega}&\qquad\text{(Parks McClellan)}\tag{7.104}\end{align}$$

where $\delta$ is the maximum allowable approximation error (Kaiser), and $\delta_1$ and $\delta_2$ are the maximum pass band and stop band errors, respectively (Parks McClellan).

The equation numbers above are taken from the book Discrete-Time Signal Processing by Oppenheim, Schafer, Buck, 2nd edition.


compute image contrast I need to get information about contrast of an image, and I compute its value in this way:

contrast=Pmax/Pmin

where Pmax is the pixel of the image with the highest value, an Pmin is the pixel of the image with the lower value.

Is this the correct way, or there is another approach?

There are many approaches to achieve this. For example, check the following work by Pederson et. al.:

Evaluation of contrast measures in relation to observers perceived contrast Marius Pedersen1, Alessandro Rizzi2, Jon Yngve Hardeberg1 and Gabriele Simone1.


will always periodogram produce correct frequencies? let say we have  some data,some signal written in this form

$x(t)=s(t)+\epsilon(t)$

where  $s(t)$ is pure signal and  $\epsilon(t)$ is  white noise  with mean zero and some variance  ,question is that if we have enough sample data,should  periodogram always  estimate frequencies  correctly?let say we have following sinusoidal data

$y[t]=A_1(sin(\omega_1*t+\phi_1)+A_2*sin(\omega_2*t+\phi_2)+....+A_p*sin(\omega_p*t+\phi_p)$+$z(t)$

i know that there is problem of variance increasing related to periodogram,but will it affect on frequency estimation?also i know that  resolution of periodogram is $1/N$,but let us suppose that  two frequency are not so close to each other,then will periodogram produce correct frequency estimation?variance is not important for me know

A periodogram created from a finite number of samples using an DFT is not really a "correct" frequency estimator, as much as it is a (leaky) frequency energy binning sieve, where the width of each bin is set by the DFT length and the window function applied.

The width and spacing of each bin may or may not be narrow enough for you to consider an error of up to half that size to still result in a useful frequency estimate.  Various interpolation methods that post-process DFT results may allow more precise estimates of spectral peaks than the bin spacing.


Discrete Time Processing of Continuous Time Signals LTI Concept I wonder if the input signal (CT) violates Shannon-Nyquist Theorem for a given sampling rate, is there any chance for the overall system not to be LTI although discrete time system is LTI?

Thanks.

Assuming a system like this:

x(t) -> C/D -> DSP -> D/C -> y(t)

where $x(t)$ and $y(t)$ are sufficiently bandlimited and the C/D and D/C blocks are ideal, then the system behaves as if it were LTI. Note that the system is not truly LTI, since it is not LTI for all inputs.

In particular, if the process of discretizing $x(t)$ produces aliases, then the output $y(t)$ can potentially contain frequencies not present in the input. Since an LTI system cannot create new frequencies, then in conclusion the system is not LTI.


How do I implement the CDF wavelet inverse? I understand and have implemented a 2D discrete wavelet transform that uses a Haar wavelet to process images.

Now I would like to implement a 2D Cohen-Daubechies-Feauveau wavelet transform, but I am uncertain about how to do the inverse transform.

The Haar transform uses a 2x2 kernel, so I have 4 coefficients for each 2x2 block of pixels, and it's obvious to me how to reverse transform to reconstruct my pixels.

The CDF wavelet uses odd tap filters, so at each image location (i,j), I have 4 coefficients.  It seems, therefore, that there is no need to downsample -- I can simply get 4 coefficient images, where each is the same size as the original image, and then run the reverse transform on those four images to reconstruct the original.  Is that correct?

This is the method I tried, but it does not restore the image -- it just gives me garbage.

If A is the low pass filter and B is the high pass filter, x is the image, and x(A,B) means apply filter A to x horizontally and B to x vertically:

x1 = x(A,A)
x2 = x(A,B)
x3 = x(B,A)
x4 = x(B,B)
y1 = x1(B,B)
y2 = x2(B,A)
y3 = x3(A,B)
y4 = x4(A,A)

Sum y1, y2, y3, y4 to get original image, x.

Does this method look correct?

I am certain that I am doing the forward transform correctly, but I am unsure about the reverse transform.  I have read the Wikipedia article on Fast Wavelet Transform, and I believe that I am implementing what this page states (without the downsampling).

To the contrary, the odd tap filters allow you to keep everything perfectly symmetric. If you make the filter centered at $0$, then the application of the filter at position $i$ requires a symmetric range of coefficients around $i$.

All 2-channel wavelets will produce a quad-tree decomposition of a 2D image, regardless of tap length. The only complication of longer, non-trivial wavelet filter banks is the treatment of the boundaries. Here one has an advantage in using symmetric filters. (Symmetric and orthogonal wavelet filters exist for 3-channel, 4-channel wavelets and so on. See articles by Heller, Belogay/Wang, Turcajeva.)

If you continue the image by mirroring at the boundaries, the transformed image will also have odd or even symmetry at the boundaries so that the size of the transformed image can be kept equal to the size of the original image.

In the reconstruction process, you have to reconstruct those mirrored coefficients in the bands and then apply the centered filters to obtain the original image back.


Discrete Fourier transform of a finite length signal which saturates at non-zero value I am performing spectral analysis of a finite length signal that saturates to a non-zero value. The signal ($s(t)$) can, practically, be write as
$s(t) = f(t) \big(1-H(t-t_0)\big) $, where $t_0$ is the time at which the signal ends and $H(t)$ is the Heaviside step function. The sudden step from a non-zero value to zero at $t = t_0$ induces artifacts in the analysis.
I am looking for advice/literature on how to deal with signals of this type.
Question: How does one deal with artifacts cause by signals terminating at a non-zero value when
performing spectral analysis.

You can & should window your signal. Choice of window depends on the specific requirements of your application. Wikipedia has a really good overview on properties and trade-offs: https://en.wikipedia.org/wiki/Window_function
Keep in mind that you already ARE windowing: your finite length window is simply an infinite length signal multiplied with a rectangular window.


Using oversampling to increase resolution of a DC-signal as input Currently I'm working on a project which uses oversampling to increase the resolution of a 12 bit ADC to a maximum of 16 bits. My goal is to fully understand the theory behind oversampling and why it is increasing the resolution. As far as I understood this topic, oversampling and decimation increases the resolution because the white noise of the input signal is distributed along a larger frequency span. After the oversampling the signal gets low pass filtered, so that we achieve fewer noise in the frequency span of interest (see the first figure).

So far so good, but I think i'm still not 100% confident of what is going on here:

How is the low pass filtering achieved? I know that I need to sample the input signal $4^n$ ($n$ = additional bits) times for every additional bit. After adding all samples together, the sum gets right shifted by $n$, which is equal to divide the sum by $2^n$. Is this the low pass filtering or do i need to somehow low pass filter it after i shifted the sum? In the application note AN118 that I used to understand this topic it seemed like the process of shifting is the low pass filter.

My input signal is a dc signal (sensor output) and I'm getting kind of a headache to understand why oversampling can increase the resolution of dc signals. Inside the appendix A of the application AN118 is a good approach to explain the reason for the increased resolution with oversampling. It also shows where the often used equation f_os = 4^n * f_s (f_os = oversampling frequency; f_s = sampling frequency) comes from. My problem is this whole appendix refers to AC-Signals... You can see this with for example equation 8 inside the application note, where the in-band noise power is calculated with the integral from 0 to f_m (f_m = highest spectral component of the input signal). In case of a dc signal as input i got f_m = 0. Am I missing something or can someone maybe explain to me why the resolution of a dc signal can be improved with oversampling?

Is something like SNR a good specification to describe the quality of a dc signal? I've read from a source that SNR is an AC specific attribute, for example in Table 2 of this Texas Instruments document. This document also states that ENOB can't be calculated for dc signals and that i need to calculate the "effective resolution" for dc signals. This was the first time i've read something about this... Can someone verify that for sure dc signals can't be specified with SNR and ENOB?

At the end of my project I would like to compare different resolutions achieved with oversampling. What would be a good specification to compare the measured results? Maybe the "effective resolution" or the overall noise calculated with the variance?
I hope someone can answere my questions. At the moment i'm struggling a bit to fully understand this topic...

My input signal is a dc signal (sensor output) and im getting kind of a headache to understand why oversampling can increase the resolution of dc signals.

The ADC puts out integers.  So, let $x$ be the integer that would come out of the ADC (100.3, say, or -333.3).  Now let $y = \lfloor r \rfloor$ be the quantization operation, where you get straight integers, like 100 or -334.
Now, let $r$ be the sum of some actual value $x$ and some nice well-behaved zero-mean noise -- $r = x + n$.  For the right sort of noise (zero-mean Gaussian that's big enough is one kind) the ADC output is a random variable with a mean equal to $x$.  So treat $y = x + n_{ADC}$, where $n_{ADC}$ is a zero-mean noise process.
Now when you average a bunch of samples of $y$, the effect of $n_{ADC}$ on the signal diminishes (in fact, its deviation goes down as $\sqrt n$ for an $n$-sample average).  That's why oversampling and averaging reduces the noise and increases resolution.

How is the low pass filtering achieved?

In the method you outline, you're averaging the input samples.  That's effectively low-pass filtering.
Note that the app note gets it slightly wrong: it's having you add up $2^n$ samples, then shift by $n$, effectively dividing by $2^n$.  However, by virtue of the averaging, you're effectively adding bits to the ADC resolution*.  Usually when I do this I just take the sum of the $2^n$ samples, knowing that there's a lot of noise in the LSBs.  If I do shift down, I shift down by $n/2$, not by $n$.

Is something like SNR a good specification to describe the quality of a dc signal?

Not really.  SNR is good when you have a well-defined signal strength.
For DC, it's much better to use $\pm x_e$, where $x_e$ is the anticipated error bars, or effective number of bits, or something.  There's plenty of examples out there -- I'd take my guidance from voltmeter specs, or from industrial instrument specs.
* You are adding to precision, not accuracy.  Accuracy is a different thing**.  The accuracy of an ADC is limited by things other than its noise characteristics.
** The best short description I have is throwing darts: accuracy is how well-centered the darts are around the bullseye -- precision is how tightly clumped they are.  If all your darts land on a point 3" over and 1" down from the bullseye, you're very precise -- but not terribly accurate.


Bandpass filtering an audio signal with MATLAB's $\tt kaiserord$: unknown new frequency component I am trying to filter an audio signal with kaiserord in MATLAB, but I'm having a problem. Here is my code:

fcuts=[50 100 20E+3 20.1E+3];
mags=[0 1 0];
devs=[0.05 0.01 0.05];
[n,Wn,beta,ftype]=kaiserord(fcuts,mags,devs, fs)  %fs from my audio signal
n=n+rem(2);
hh=fir1(n,Wn,ftype,kaiser(n+1, beta),'scale');
filtered=filter(hh,1,W) %W is my audio signal

I try to filter audio with $50\textrm{ Hz}-20\textrm{ kHz}$ bandpass filter but after applying filter; I realize that there is signal before $47\textrm{ Hz}$,

Where is my mistake?

Your filter is a stop band attenuation of about 40 dB. Frequencies outside your passband will be attenuated by at least 40 dB but they will not go down all the way to zero.

The FIR filter is very long (2000 or samples). It takes that many samples for the filter to reach "steady state". In other words, you should throw away the first 2000 samples of your filtered signal before doing any analysis on it.


general n-dimensional name of Pixel and Voxel In 2D images we call the smallest picture elements pixels, in 3D images we call them Voxels (for volume pixels).

I was wondering if there is a word that can denote them both (pixel and voxel), and possibly be used to mean any n-dimensional picture element.

For example, I want to write something like:

  Given an n-dimensional image with N pixels/voxels

Where I want to replace pixels/voxels with just one word that would also denote these elements in other dimensions.

From Wiktionary:

  hypervoxel (plural hypervoxels)

  A multidimensional analogue of a voxel

From Digital Design Media - Mitchell & McCullough, 1995:


signal of $\frac{1}{f}$ Noise I am working on kind of signals that can be classified as 1/f Noise.
In this article is explained that the time history of 1/f noise follows a power law proportional to $t^{-\alpha}$. I am very surprised about this, because when I think to this kind of signal histories, they do not seems to me proportional to $t^{-\alpha}$ or a series of $t^{-\alpha}$.
Have someone an explanation for this?
I have also read this article, that highlights the proportionality to $t^{-\alpha}$ and introduces the theory of the "Self-Organized criticality". But I didn't find an explanation to my problem.

Power-law behaviors in frequency can be found in several unrelated observations and systems. This is apparently the case for $1/f$ or flicker noise. Note that an exact $\alpha=1$ exponent might be too stringent, and people can be interested in a wider range, like  $1/f^\alpha$ with $1/2\le\alpha\le3/2$.

The plot above is from 1/f noise.
Since only looking at the decay of a PSD does not provide you with much ideas about the morphology of signals, you might be interested in looking at generative processes, or accumulated actions that may generate such types of spectra, from physics phenomena (superposition of relaxation processes, infinitely large fluctuations, noise in diffusion processes etc.). And then maybe find analogies when the phenomenon you are studying.
Here are three papers of interest:

$1/f$  noise: a pedagogical review by Edoardo Milotti
$1/f$ noise’’ in music: Music from $1/f$ noise by Richard F. Voss and John Clarke

The spectral density of fluctuations in the audio power of many
musical selections and of English speech varies approximately as $1/f$
(f is the frequency) down to a frequency of 5×10−4 Hz. This result
implies that the audio‐power fluctuations are correlated over all
times in the same manner as ’’$1/f$ noise’’ in electronic components.
The frequency fluctuations of music also have a $1/f$ spectral density
at frequencies down to the inverse of the length of the piece of
music. The frequency fluctuations of English speech have a quite
different behavior, with a single characteristic time of about 0.1 s,
the average length of a syllable. The observations on music suggest
that $1/f$ noise is a good choice for stochastic composition.
Compositions in which the frequency and duration of each note were
determined by $1/f$ noise sources sounded pleasing. Those generated by
white-noise sources sounded too random, while those generated by
$1/f^2$ noise sounded too correlated.

Power laws, Pareto distributions and Zipf's law


How to remove ECG artifacts from EMG data? I'm new to signal processing and try to analyze EMG data from pectoralis major and other respiratory muscles. Unfortunately, there are strong ECG artifacts in all recordings and I don't know how to remove them from my data.

My data looks like this:

In my experiment, there are two expected EMG-events, that should be at around second 2 and 8 of the recorded data:

I'm working with python and matlab and couldn't find code-examples for this specific problem.

My sampling rate is 4000 Hz. Here is an example of the EMG-data.

Thanks a lot for your help!

Here is principal component analysis of the data you supplied. PCA will extract the features within the signal as is done in picture 1.

Note that the signal will be broken down as it is presented in picture 2 from the paper in 3 (page 15). To remove, just note that most of the principal components look like they come from the ecg data, so lets just take all the eigvectors after 10 and project them onto our signal and add them up to get emg data (see picture 4). I'm not entirely sure what the emg data is suppose to look like, but this looks close to diaphragm emg data I see on google images. Try excluding even more eigenvectors to see if you get what you are looking for.

    close all
    clear all
    H = csvread('\example.csv');

%regroup data into snapshots of the signal
%to treat as a group of random processes
l = 2710

Z = (H(1:l*10));
X = reshape(Z,[l, 10]);[![enter image description here][5]][5]
RV = X*X'/length(X)
%eigen decomposition
[U,D,V] = svd(RV);

figure
stem(diag(D))
title('eigenvalues')

%projection
figure
for i  = 1:12
xx0 = V(:,i).*Z(1:l);
subplot(4,3,i)
plot(-xx0)
strr = sprintf('principal component: %d', i);
title(strr)
end

sig = 0;
for i = 11:length(V)
  sig = sig + V(:,i).*Z(1:l);
end

figure
plot(sig)
title('Emg data?')

figure
subplot(121)
plot(V(:,11).*Z(1:l))
title('Principal component 11, emg?')
subplot(122)
plot(V(:,12).*Z(1:l))
title('Principal component 12, emg?')


How to plot histogram of difference of two images using MATLAB i would like to calculate difference of two images. Then, histogram of difference of two images. So I can plot a histogram but I am not sure whether it is true or not? Also, I am not sure about what should be second parameter of built-in histogramfunction.

My two images are original_result and output_image. And difference of two image is image_of_dif.
my result:

image_of_dif = double(original_result) - double(output_image);
figure; histogram(image_of_dif(:), 511);

In your second line, casting with uint8() will change a negative value to 0. Avoid casting or change to another type that supports negative values.


Question about Hilbert transform of a cosine signal In the paper entitled "On instantaneous frequency", it is claimed that the Hilbert transform of a signal in the form of $A(t)cos(\phi(t))$ would result in $A(t)sin(\phi(t))$, where $A(t)$ is instantaneous amplitude and $\phi(t)$ is the phase function.
As an example, consider the following:
$$t = [0, 20],\;\textrm{time interval in seconds, sampled each 0.01 sec}$$
$$A(t) = e^{-0.1(t - 10)^2},\;\textrm{instantaneous  amplitude}$$
$$\omega = [0.1,\ 0.5,\ 1.0,\ 1.5],\;\textrm{angular  frequency  in rad/sec}$$
$$f(t) = A(t)\cos(\omega t)$$
$$g(t) = A(t)\sin(\omega t)$$
$$h(t) = \frac{1}{\pi}P\int_{-\infty}^{+\infty}\frac{f(\tau)}{t - \tau}dt,\;\textrm{Hilbert  transform  of } f(t)\textrm{,  where } P  \textrm{ is the   Cauchy  principal  value}$$
It is expected that $g(t)$ and $h(t)$ be the same according to the fact claimed in mentioned paper. I used function "hilbert" of MATLAB to compute Hilbert transform of $f(t)$ for aforementioned $\omega$ values. Results show that for $\omega$ close to zero, $g(t)$ and $h(t)$ are different and as $\omega$ gets higher, the difference decreases. I was wondering if I am making a mistake or something is wrong with function "hilbert" of MATLAB?
The following figure shows the results:

Make sure that you understand the conditions under which

$$\mathcal{H}\big\{A(t)\cos(\omega_0 t)\big\}=A(t)\mathcal{H}\big\{\cos(\omega_0 t)\big\}=A(t)\sin(\omega_0t)\tag{1}$$

holds. Eq. $(1)$ holds if $A(t)$ is a low-pass signal with a cut-off frequency smaller than $\omega_0$. This implies that $A(t)\cos(\omega_0t)$ is a band-pass signal with no energy at DC. Eq. $(1)$ is a special case of Bedrosian's theorem.

In your example, $A(t)$ is not band-limited, and, consequently, $(1)$ doesn't hold. However, if you increase the modulation frequency $\omega_0$, the overlap of the spectra of $A(t)$ and of the carrier becomes smaller, because $A(t)$ does have a low pass character. So for large $\omega_0$, Eq. $(1)$ is approximately satisfied, which is exactly what you see in your plots.

Also take a look at this related answer.


Least Squares Fitting to Inverse Exponential Function I have a time series of measurements that resembles the shape of an exponential function. The samples are a bit noisy and sometimes there is a weak sine like ripple signal ontop of it.

Simplified the function looks like this:

$$ y = a (1-e^{-bt + c}) + d $$

The data that I have lies between $t=$0 to 2, e.g. where the function is still curvy.

My goal is to find $a+d$, e.g. the value of the asymptote or the $y$ value of the graph above at infinite time.

My idea so far is to take the samples that I have, do least square fitting on the function above and then just solve it.

Q: Will this approach work or is there a easier way?

Q: I'm a bit lost how to do a least square fitting to a function like the one above. With a falling decaying exponential I would just take the log and least square fit a line, but how do I do this with my rising function?

Oh, one more thing: I don't need to do this as a one-off so I can't just plug the data into a general purpose solver. I have to implement this later for resource constrained microcontroller.

I think you are going to be stuck with an iterative approach.

$$ y = a (1-e^{-bt + c}) + d $$

First rewrite your equation like this:

$$ y = g + h e^{-bt} $$

Where:

$$ g = a + d $$

$$ h = -a e^c $$

You have a set of values $(t_n, y_n)$.  Consider the $y$ values as a single vector $\vec Y $.  Then pick an initial value for $b$ and construct a vector $\vec B$ using your $t$ values in $ e^{-bt} $.  You now have a vector form of your problem:

$$ \vec Y  = g \vec U  + h \vec B $$

Where $\vec U$ is a vector of all ones.

You can solve for $g$ and $h$ by dotting $\vec Y$ with $\vec U$ and $\vec B$.  This gives you a linear system of two equations two unknowns which is easily solvable.

$$ (\vec Y \cdot \vec U)  = g (\vec U \cdot \vec U) + h (\vec B \cdot \vec U) $$

$$ (\vec Y \cdot \vec B)  = g (\vec U \cdot \vec B) + h (\vec B \cdot \vec B) $$

Note that $ \vec U \cdot \vec U $ is your sample count and $ \vec Y \cdot \vec U $ is the sum of your $y$ values.  Calculating $ \vec U \cdot \vec B $, $ \vec B \cdot \vec B $, and $ \vec Y \cdot \vec B $  can be done in single loop without actually having to store the vector values if you are memory constrained. Storing $ \vec B $ will save recalculation when you are measuring the fit.

To measure your fit, define the difference vector:

$$ \vec D  = \vec Y  - \left( g \vec U  + h \vec B \right) $$

Your square fit is now:

$$ f = \vec D \cdot \vec D $$

Your problem now is to find the value of $b$ which minimizes $f$.  There are several ways to do this.  The fastest converging one might be to pick three values of b, calculate an f for each, fit a parabola to your $(b,f)$ points, find the minimum $b_*$.  To iterate, choose three new tighter $b$'s around $b_*$.  Repeat until good enough.

The value of $g$ is the answer you are looking for.

Hope this helps,

Ced

Followup

As I mentioned in my comment to Ben, this equation is overparameterized.  My recommendation would be to get rid of $d$ by setting it to zero.  This leaves

$$ a = g $$

$$ c = \ln( -h/g ) $$

As a differential equation, per r b-j:

$$ y = g + h e^{-bt+c} $$

$$ y' = h e^{-bt+c} (-b) $$

$$ y' = b ( g - y ) $$

Getting a good initial value of b:

There's nothing like sleeping on a problem.  If your samples are evenly spaced in time then you can find a good initial value of b.

The simplest involves using three evenly spaced points at time $t_A$, $t_B$, and $t_C$.

$$ t_B = t_A + s $$
$$ t_C = t_B + s $$

Calculate the following ratio:

$$ r = \frac{ y_B - y_A }{ y_C - y_B } $$

The value of b can be now be solved for:

$$ r = \frac{ ( g + h e^{-bt_B} ) - ( g + h e^{-bt_A} ) }{ ( g + h e^{-bt_C} ) - ( g + h e^{-bt_B} ) } $$
$$ r = \frac{ h e^{-b( t_A + s ) } - h e^{-bt_A} }{ h e^{-b( t_B + s ) } - h e^{-bt_B}  } $$
$$ r = \frac{ e^{-b t_A }( e^{-bs} + 1 ) }{ e^{-b t_B }( e^{-bs} + 1 )  } $$
$$ r = e^{-b ( t_A - t_B ) } $$
$$ b = \frac{ \ln( r ) } { t_B - t_A } $$

In order to mitigate the effects of noise, more points can be used.  Another formula would be:

$$ r = \frac{ y_1 - y_2 + y_3 - y_4 + y_5 - y_6 }{ y_7 - y_8 + y_9 - y_{10} + y_{11} - y_{12} } $$

$$ b = \frac{ \ln( r ) } { t_7 - t_1 } $$

Any number of formulas can be created.  The rules are that the coefficients have to add up to zero and the time spacing has to be the same in the numerator and the denominator.  Also, you want to stay on the steeper part of the curve so the differences in values are larger compared to your noise level.

With enough sample values used, the value of $b$ might be good enough that the first calculation of $g$ will be accurate enough that no iterations are needed.


FFT Resolution Inference Let's say I have an audio clip of $N$ samples.  I want to describe every $10$ ms of audio with an FFT of size $1024$, without information bleeding over from other $10$ ms windows.  Given that the sample rate of this audio clip is a standard $44100$, I can't do that with a regular FFT.  A regular FFT would give me $441 / 2 = 220$ bins.

Is it possible to take many, many FFTs of differing sizes over different windows to infer a larger frequency resolution for a smaller timeframe?  Let's say I take a $1024$-FFT for each sample range increasing by $t = 1$ (samples $0$ through $2047$, then samples $1$ through $2048$ and so on), then a bunch of $512$-FFTs, then a bunch of $256$-FFTs.  Would this large dataset contain enough useful information to "subdivide" the bins of the smaller, $220$-FFT?

You can use interpolation (or equivalent zero-padding) to get your desired sub-divided 512 frequency points from your 220 FFT bins.  Use a Sinc interpolation kernel for a rectangular window, or the transform of the window for other windows.  If you want more real frequency resolution (as in clear separation of multiple peaks present in the frequency spectrum), you will have to use (overlap or "smear") more than 10 mS of data (perhaps over 50 mS).


Smoothing in frequency domain and impulse reponse of the filter in time-domain I have a low-pass signal in frequency domain, $$\left[f_{min}=0, f_{max}=20\,\mathrm{GHz}\right]_{n=100}$$

After re-configuring the input signal as
$$X=\left[-f_{max},\ \dots\ -f_{min},\ \mathrm{DC},\ f_{min},\ \dots\ f_{max} \right]_{2n+1}$$

and constructing a Hamming window of length $2n+1$ as win = hamming(2n+1),

Now the question is, how to perform smoothing and eventually get the impulse response of the filter?

Another confusing point is whether I should convolve the input signal with the window function?

-or-

Is it the multiplication between the input signal and the window function? Like X.*win?

In this context do I need perform fft(hamming(length(2n+1))) on window function?

time            = 5; % in nano-seconds
z0              = 50;
A               = importdata("sparameters_lowpass.dat");
M               = 2048;

freq            = A(:,1);
freq_neg        = -1.0*flip(freq);

realPart        = A(:,2);
imagPart        = A(:,3);

s11_freq        = realPart + imagPart*i;
s11_conj        = conj(flip(s11_freq));

F               = [freq_neg(1:end-1)', freq(1),     freq(2:end)'];
X               = [s11_conj(1:end-1)', s11_freq(1), s11_freq(2:end)'];

win             = (hamming(length(X)));

freq_multip     = X.*win';

s11_time        = ifft(freq_multip);

z_in            = z0 * (1 + s11_time)./(1 - s11_time); % input impedance

Any help to the correct approach is very helpful.

You get the time domain values of hamming window from the hamming function. Multiplying it with the frequency domain data will not give the correct result. Smoothing is usually performed by multiplying in time domain. You could try multiplying the fourier transform of hamming window with the Frequency domain signal. and then performing the inverse fft.


How do I calculate the eigenvalues of a covariance matrix which contains harmonic functions I have read that $ce^{-\frac{\phi_i-\phi_j}{\rho}}$ is a harmonic function of the form $e^{-in\phi_i}$, and therefore it's eigenvalues are one of the Fourier components of $(|\phi_i-\phi_j|)$.

My questions are:

How come $ce^{-\frac{\phi_i-\phi_j}{\rho}}$ has the form of $e^{-in\phi_i}$ if it is a real function?
How do I calculate the eigenvalues of such a function, or, if it is a difficult calculation, is there a known theorem(s) I can base this result on?

I am not an expert in harmonic analysis, but I would like to know more about this subject (without getting into too much details). I need it to my research in neurosciences.

Here is a link to the article I am reffering to: Population coding

Thanks!

The included document does not say $ce^{-\frac{\phi_i-\phi_j}{\rho}}$ has the form of $e^{-in\phi_i}$. The included document uses $C$ for three completely different things:

$C_{ij}$ is the correlations between $i$ and $j$ items.
$C(\phi_i-\phi_j)$ a component of that correlation, dependent on the angles between the two items $i$ and $j$.
$C_n$ the eigenvalues of the covariance matrix made up of the $C_{ij}$.

The term $e^{-in\phi_i}$ is a component of $C_n$ the eigenvalues / Fourier components expression in (31).

The idea would be to form the matrix $C_{ij}$ and then take the eigensystem decomposition of that matrix.


convolve to differentiate black and white colors a figure for instance of size 500*500 has half above part with black and below half white should result in a white line where the
white meets the black (something like  a single line at line 250 with pixels as 255) on convolution with a matrix.
 I am trying to make a matrix, but not able to visualise how the matrix would be? can some one help me build the matrix for this

You have to make all the values zero in the matrix except the middle row and thus, what you can start thinking is how to get such filter matrix which can give you this result. One of the many filter matrices you can think of, which satisfy your requirements is [0 1 0;0 0 0;0 -1 0].


Is there a known correlation between human voice (360Hz to 5kHz) and the frequencies above/below that? Hello SignalProcessing!

Given the important frequencies of a human voice (360Hz to 5kHz) to communicate and the frequencies above and below that make the voice sound way better.

Imagine we are recording a phone call... the provider normally transmitts only the frequencies needed for moderate communication. I am interested in reconstructing the information for the lost frequencies above and below by training deep neuronal nets.
To improve the training progress it would be helpful to figure out if there is sort of a correlation.

Any reference about that topic would be great!
Thank you.

There is a correlation and it can be learned with machine learning methods. A simple google search can give you many papers, recent research

Waveform Modeling and Generation Using Hierarchical Recurrent Neural Networks for Speech Bandwidth Extension
Zhen-Hua Ling, Yang Ai, Yu Gu, Li-Rong Dai

older paper:

Techniques for the Regeneration of Wideband Speech from Narrowband Speech
Jason A. Fuemmeler, Russell C. Hardie, William R. Gardner


Is the value of a discrete time signal 0 when the index is a non integer value? Is it safe to assume that if the index of discrete signal is not an integer, then the value of it is 0? For example is the following true:

If we have $x[n] = 1 + (−1)^n$ does $x[3/4]=0$?

The TA in my Digital Signals class is saying that this notion is correct, however the general consensus from what I've seen from Google is that $x[3/4]$ should be undefined, not zero.

Discrete-time signals are, unsurprisingly, only defined for discrete times, i.e. for integer values of their argument. The value of $x[n]$ for $n$ not an integer is undefined, and $x[3/4]$ or $x[\pi]$ is meaningless.


Finding the Best Gaussian Smoothing Kernel to Minimize the Discrepancy Between Two Images Suppose we have two grayscale images, $A$ and $B$. $A$ and $B$ very strongly resemble each other, such that the mean of the absolute difference $\lvert A - B\rvert$ is fairly small. Suppose further that $B$ actually appears to be a blurred version of $A$, although it was not necessarily obtained through a blurring process. I'd like to find the ideal Gaussian kernel $g$ to smooth $A$ with to maximize its correspondence to $B$, in other words, select $g$ s.t.

$$\underset{g}{\arg\min}\sum{\lvert(g*A) - B\rvert}$$

where $g$, $A$, and $B$ are 2D arrays and $*$ is the convolution operator.

This is closely related to Blind Deconvolution.
The only difference is we limit our self to a very specific type of blur kernels.

The nice thing about the Gaussian Kernel is being defined by single parameter - The Standard Deviation of the kernel.
The less nice thing is the connection isn't linear.

Optimization Problem

Let's define a classic non linear model for this problem:

$$ \arg \min_{\sigma} \frac{1}{2} {\left\| A \left( \sigma \right) x - b \right\|}_{2}^{2} $$

Where $ A \left( \sigma \right) $ is the convolution matrix generated by a Gaussian Kernel parameterized by $ \sigma $, $ x $ is the original image ($ A $ in your question) in a vector shape, and $ b $ is the blurred image ($ B $ in your question).
This is a classical Non Linear Least Squares problem which can be solved by MATLAB using lsqnonlin().

Code Sample

This is the main part of the code:

mA = rand([numRows, numCols]);

gaussianKernelStd = kernelStdLowerBound + ((kernelStdUpperBound - kernelStdLowerBound) * rand(1));
mB = ApplyGaussianBlur(mA, gaussianKernelStd, STD_TO_RADIUS_FACTOR);

% Objective Functions
hObjFun = @(kernelStd) reshape(ApplyGaussianBlur(mA, kernelStd, STD_TO_RADIUS_FACTOR) - mB, [numPx, 1]);
estKernelEst    = lsqnonlin(hObjFun, initKernelStd, kernelStdLowerBound, kernelStdUpperBound, sSolverOptions);

Results

As can be seen, the estimation is almost perfect.
On large images it might take time (Using some tricks of the Gaussian Filter in the Fourier Domain one could do that there with major speed up), but still, it is not free.

The full code is available on my StackExchange Signal Processing Q49121 GitHub Repository.


How do we impose Hermitian Symmetry in optical OFDM? I’m currently trying to wrap my head around optical OFDM. I’m reading various tutorials and articles, but I’m hitting a stumbling block trying to understand Hermitian Symmetry and how it is applied.
I’m taking the block diagram of an OFDM system from Figure 3 in this paper, and working through it step by step starting from the input to the transmitter. It’s probably best that I go step-by-step up until the point where I get confused.

Step 1 - serial to parallel
For the sake of simplicity I’m skipping the “Coding” and “Interleaving” steps from that paper. Let’s start with a serial bit stream we want to send: 1, 1, -1, -1, 1, 1, 1, -1, 1, -1, -1, -1, -1, 1, -1, -1, -1, 1…etc.
We will map this to parallel bit streams, one for each subcarrier. I’m using 4 subcarriers in this example, C1 to C4. The parallel streams are now:

C1
C2
C3
C4

1
1
-1
-1

1
1
1
-1

1
-1
-1
-1

-1
1
-1
-1

-1
1
1
-1

-1
-1
1
1

Step 2 - QAM modulation
Modulate each subcarrier. Here I'll use 4-QAM. Taking 2 bits at a time from each of the subcarrier columns I get:

C1
C2
C3
C4

1 +1j
1 +1j
-1 +1j
-1 -1j

1 -1j
-1 +1j
-1 -1j
-1-1j

-1 -1j
1 -1j
1 +1j
-1 +1j

Step 3 - IFFT
This is where my understanding starts to break down. I understand for “normal” OFDM we’d pass these values to the IFFT to generate the time-domain OFDM signal. However, what exactly are we sending to the IFFT block? Do we send it each row of that table in turn, i.e. the values of each subcarrier at a particular sample time?
Step 4 - Hermitian Symmetry
This is where my understanding really breaks down! I understand that Optical-OFDM require purely real-valued time-domain OFDM signals. I also follow that if the input to the IFFT block has the property of "Hermitian symmetry", then the IFFT output will be the desired real signal.
Hermitian symmetry is "imposed" on Xk symbols by ensuring that:
X0 = XN/2 = 0 and
Xk = X*N-k for 0 < k < N/2
Where N is the IFFT length and k is the carrier number (I assume?)
My confusion here is if we have to "impose" that certain carriers must take a given value (e.g. carrier 0  and carrier N/2 must equal 0) then isn't that somehow "overwriting" the value for that carrier that we have applied at the QAM modulation stage? That seems like it must be wrong to me, but I can't understand why!
So, apologies for the long-winded question, but just to clarify I'm asking if my understanding of the IFFT at Step 3 is correct, and what exactly we are doing when we "impose" Hermitian Symmetry onto our signal prior to the IFFT in Step 4.

I don't understand your question about the IFFT step: as in "normal" OFDM, you send a full vector to the IFFT, as that is a vector-valued operation of vectors. The output are (interpreted as) time-domain samples.
In fact, you really shouldn't have put Step 3 and 4 as separate operations, and certainly not in that order: the input to the IFFT has to be hermitian symmetric to yield a real-valued signal (which is what you need if you don't want to implement a quadrature upconverter to mix your signal into a carrier frequency, so in all honesty, fig. 3 is misleading in that the $-\sin(2\pi f_c t)$ branch doesn't have anything to do if you ensure your IFFT input is hermitian symmetric.)

My confusion here is if we have to "impose" that certain carriers must take a given value (e.g. carrier 0 and carrier N/2 must equal 0) then isn't that somehow "overwriting" the value for that carrier that we have applied at the QAM modulation stage?

No, you're just mixing up what you're doing: you're not free to assign an arbitrary value to each subcarrier at will, you need to make sure the vector is hermitian symmetric. You hence only have half the subcarriers to choose freely, the other one is inherently set to the complex conjugate.
So, you're not "overwriting" anything. You're just only assigning your QAM symbols to half the subcarriers.


Compute Similarity between Fourier Transforms I'm looking to compare the Fourier Transforms generated by accelerators and gyroscopes that collected data of people walking. I've looked to see if there is a standard form of comparison, but I have yet to find one.

Here's my specific use case: I'm looking to compare similarity between the dominant signals of two Fourier Transforms (i.e. how close are they to being the same frequency/magnitude). Is there a specific metric that would help identify that?

If not (or perhaps irrespective of my goal), what are some common ways of comparing two Fourier transforms? I've found surprisingly little on the subject.

Most simply, you can compute simple statistical measures on each FFT frame and compare those directly. The signal energy, its distribution, variance, etc across the frequency spectrum are interesting measures. You can also directly compute the correlation coefficient between accelerometer signal pairs (normalized to their length).

It might be useful to look into research on “gait recognition”, the literature points to “low frequency-domain entropy” as feature that is useful in discriminating accel/gyro data that differ in “complexity”

Have you seen this thread on comparing 3d data from accelerometers. It is worth considering if using the raw 3d data, doing a dimensionality reduction (PCA, etc) and then using a simple distance metric may get you further than analyzing raw FFT output. Doing an FFT as you've proposed doing is only one way to simplify the data, but I think the Fourier transform is too coarse & noisy to give you directly comparable results, every time.

As an alternative, it might be easier/faster to use machine learning, i.e. training a "people walking" classifier using any modern supervised learning technique.


Estimate Parameters of Linear Combination of Harmonic Signals with Partial Known Frequencies Consider real signal of the form:
$$
\begin{align*} x(n) = &\sum_{m=1}^{M} \left( c_m \cos(2\pi f_m n) + s_m \sin(2\pi f_m n) \right) =\\
&\sum_{m=1}^{M} a_m \cos(2\pi f_m n + \phi_m)
\end{align*}
$$
where signal values of $x(n)$ are known for $n = 1,\dots, N >> M$ and some frequencies $f_m$ are known real numbers in $(0.0,0.5)$.
What is the best way to estimate $c_m$ and $s_m$ real parameters?
Estimation of the form:
$$c_m + j s_m = 2 \frac{\displaystyle\sum_{n=1}^{N} x(n) \exp(2 \pi j f_m n)}{N}$$
is not very accurate ($j^2=-1$), but windowing (e.g. cosine window) improves estimation accuracy significantly:
$$c_m + j s_m = 2 \frac{\displaystyle\sum_{n=1}^{N} w(n) x(n) \exp(2 \pi j f_m n)}{\displaystyle\sum_{n=1}^{N} w(n)}$$
Why window is so effective here? Is it some kind of well-known formula?
The least-squares method also gives good results if frequencies are known for all terms.
Are there any other good methods for the case when not all frequencies are known? Perhaps based on SVD or some kind of orthogonalization.
Test code in Mathematica:
F1 = 0.123000 ;
F2 = 0.312874 ;

LENGTH = 1024 ;
RANGE = Range[LENGTH] ;

SIGNAL = 0.0 ;
SIGNAL = SIGNAL + 1.0*Cos[2*Pi*F1*RANGE] + 0.5*Sin[2*Pi*F1*RANGE] ;
SIGNAL = SIGNAL + 0.05*Cos[2*Pi*F2*RANGE] + 0.01*Sin[2*Pi*F2*RANGE] ;

WINDOW = HannWindow[N[Rescale[Range[LENGTH],{1,LENGTH},{-1,1}]]]^2 ;

(* NO WINDOW *)
2*(Total[SIGNAL*Exp[2*Pi*I*F1*RANGE]]/LENGTH)
2*(Total[SIGNAL*Exp[2*Pi*I*F2*RANGE]]/LENGTH)
(* OUTPUT: 0.9994929536863852`\[VeryThinSpace]+0.49994639046951617` I *)
(* OUTPUT: 0.04819835064062154`\[VeryThinSpace]+0.012127390241918076` I *)

(* WITH WINDOW *)
2*(Total[SIGNAL*WINDOW*Exp[2*Pi*I*F1*RANGE]]/Total[WINDOW])
2*(Total[SIGNAL*WINDOW*Exp[2*Pi*I*F2*RANGE]]/Total[WINDOW])
(* OUTPUT: 0.9999999999775656`\[VeryThinSpace]+0.5000000000025541` I *)
(* OUTPUT: 0.05000000004559643`\[VeryThinSpace]+0.009999999958900167` I *)

(* LEAST SQUARES *)
(* ONE KNOWN *)
LeastSquares[{Cos[2*Pi*F1*RANGE],Sin[2*Pi*F1*RANGE]}//Transpose,SIGNAL]
(* ALL KNOWN *)
LeastSquares[{Cos[2*Pi*F1*RANGE],Sin[2*Pi*F1*RANGE],Cos[2*Pi*F2*RANGE],Sin[2*Pi*F2*RANGE]}//Transpose,SIGNAL]
(* OUTPUT: {0.999957445314386`,0.4999499273288035`} *)
(* OUTPUT: {1.`,0.49999999999999983`,0.05000000000000022`,0.010000000000000056`} *)

Bottom Line: The OP is experiencing the effects of spectral leakage and the primary motivation for windowing with the Discrete Fourier Transform. The reason the windowing approach achieves much better estimation accuracy is explained as follows:
The OP's equation is a sum of real sinusoids, and his estimator to determine those sinusoids is a form of the Discrete Fourier Transform (DFT). Without additional windowing (which is a rectangular window) it is well understood that the DFT suffers from a considerable amount of spectral leakage, in that the Kernel of the rectangular window (the Discrete-Time Fourier Transform or DTFT of the rectangular window which is of the form of an aliased Sinc function) convolves in frequency with the underlying two impulses in frequency of each real sinusoid. Since the aliased Sinc function decreases relatively gradually with frequency offset from the location of the impulses representing the tones, and has such high side-lobes, many other frequency bins are affected causing errors in their own measurements (and a significant decrease in dynamic range or the ability to detect a lower level tone in the presence of a strong one). Windows are selected that have the property of much lower side-lobe levels, and roll-off much faster in frequency, which therefore causes any given tone to have much less impact on the estimation of other tones present.
This is made much clearer if we observe the DTFT instead of the DFT, as the DTFT is a continuous function of frequency and the DFT is samples on this continuous function. By observing the DTFT we can more easily get the full intuitive picture of what is being described above. Further we also see how approximating the DTFT (through zero-padding prior to doing the DFT) can also be used to approximate the frequencies). I show this with plots of the OP's original frequencies where we can see the effect of the windowing directly:
First consider the OP's example of 1024 samples with just F1 alone with and without additional windowing. First without additional windowing (rectangular window) The log magnitude plot of the DTFT with the corresponding 1024 DFT samples of just F1 is given in the plot below:
This is specifically the DTFT of
$$sig1 = 1.0\cos(2\pi f_1 n) + 0.5\sin(2\pi f_1 n)$$
for n = 0 to 1023.  And f1 = 0.123
(Note: the above is just a single tone at $f1$ with additional phase shift as can be derived from the relationship $a\cos(x)+b\cos(x) = \sqrt{a^2+b^2}\cos(x-\tan^{-1}(b/a))$
Magnitude Spectrum of F1 Signal Only - Rectangular Window

If we zoom in on one of the spectral peaks, we can clearly see the Sinc shaped Kernel (in this case given the number of samples the aliasing would be very small so this is a very close approximation of an actual Sinc) that has convolved with the underlying impulse that would be at the F1 frequency location, and we see the spectral leakage and how it affects the other frequency locations. F1 as 0.123 is very close to bin number 126 as given by $1024 \times .123 = 125.952$.
Zoom In on F1 Peak Location

Below I show the two frequency components as the OP has used, but still as separate signals to demonstrate how the spectral leakage of each will effect the measurement of the other.
This is specifically
$$sig1 = 1.0\cos(2\pi f_1 n) + 0.5\sin(2\pi f_1 n)$$
$$sig2 = 0.05\cos(2\pi f_2 n) + 0.01\sin(2\pi f_2 n)$$
For n = 0 to 1023 and f1 = 0.123 and f2 = 0.312874
Magnitude Spectrum of F1 and F2 Signals (Separate) - Rectangular Window

Notice if we zoom in on the peak for F2 how significant the leakage from F1 is will effect the estimate of the magnitude for F2, given that the leakage from F1 is only 29 dB lower! If the magnitude of F2 was made lower, the resulting error would be even more pronounced. (Also notice how with the DTFT we can locate the F2 frequency precisely at its fractional bin location of $1024 \times 0.312874 \approx 320.3830$ which is where the peak of the orange curve in the plot below would be located, but once combined with the F2 signal given by the blue curve, this would cause a frequency error as well since the peak of the blue curve happens slightly to the right of the orange F1 curve.
Zoom In on F2 Peak Location - With Rectangular Window

If we then compute the DTFT of the composite waveform with both F1 and F2, we see the resulting amplitude and frequency error in vicinity of the F2 Peak location as shown in the plot below (where Green is the DTFT of the composite waveform summing both F1 and F2 components).

This error is predicted by the coherent sum of the two waveforms and thus can be predicted from rectangular window Kernel and knowledge of the relative amplitude and frequency offsets of the two signals. Given the waveforms are complex the phase of the components will impact the summation result (thus the inteference could feasibly reduce or increase the amplitude of the affected component).
In this example the following were the signal contributions at the F2 peak location (not in dB scale as above to show the summation directly).
F2 alone at bin 320.38: magnitude = 52.10 (given as 52.00 -j3.26)
F1 alone at bin 320.38: magnitude = 1.70 (given as 1.42-j0.94)
F1 + F2 at bin 320.38: magnitude = 53.58 (given as 53.42-j4.19)
Where F1 + F2 as measured is predicted by the sum of F1 and F2 alone when the complex components are added.
Finally if we observe the above plot with the use of windowing (we'll use the same Hann window the OP used) we can see the significant reduction in spectral leakage and then exactly how this would contributed to a better estimation of both the frequency and magnitude of the underlying signals.
Magnitude Spectrum of F1 and F2 Signals (Separate) - Hann Window

Zooming in on the F2 spectral peak now shows that the leakage from F1 is now over 100 dB lower, but we also see the drawback of windowing- the main lobe of the Kernel is wider leading to reduced frequency selectivity. (The rectangular window has the best selectivity with the narrowest possible main lobe, but the worst dynamic range due to the high sidelobe level and slow roll-off).
Zoom In on F2 Peak Location - With Hann Window

For more info on DTFT vs DFT and approximating the DTFT with zero padding see Getting the DTFT from the DFT samples and For 2D signals can it be said that the frequency response is the same as the Fourier transform?
For additional details on windowing see:
Why would one use a Hann or Bartlett window?
For additional details on spectral leakage see: Intuition for sidelobes in FFT
Suggestions for Efficient Processing
Assuming the objective is to only find the magnitudes for each of the tones, if 0.01 dB accuracy is sufficient, a suggestion is to use a flattop window, which provides nearly 100 dB of sidelobe rejection while maintaining a very flat main lobe in frequency extending out to $\pm 0.5$ bins. This means no further interpolation would be necessary to estimate the magnitude of a tone located at a fractional bin frequency to get an estimate within 0.01 dB of the actual magnitude.
Below shows an example of the Kernel (DTFT) of a flattop window with just 128 samples, showing nearly 100 dB of sidelobe rejection reached at a 5 bin offset.

Zooming in to see the flatness, we see that with 128 points we maintain less than 0.01 dB flatness over $\pm 0.5$ bins.

If the $0.01$ dB magnitude accuracy was acceptable, the solution to minimize processing would be to choose the maximum DFT size required to ensure no two tones are spaced closer than 5 bins based on the minimum tone proximimity (a 58.7 dB leakage which occurs at approximately 4.6 bins would cause a maximum amplitude error of 0.01 dB). This would set the minimum number of samples to process, after which can be multiplied by the flattop window and fft'd where the maximum bins would be sufficient to provide the amplitude estimates.
Here is an interesting paper with further information on designing the flattop window to any passband error specification (at the expense of sidelobe performance), as well as other optimized flattop solutions: Design and Improvement of Flattop Windows with Semi Infinite Optimization.


solving recurrence relation with Z-transform (LTI and causal) I just wanted to doublecheck answers for my sanity's sake (exam next week)

problem statement

recurrence relation, solve it

$y[n+1]= 35 + y[n]*0.5$

according to my teacher it will be such that the input sequence must be causal, so that we have re-arranged form such as

$y[n+1]-  y[n]*0.5= 35u[n]$
,where $y[0]=20$, and  causal input

My understanding is that we should go as follows

get $X(z)$
find what is $Y(z)$
then we can maybe do inverse Z-transform to get the y[k] output sequence in explicit form sequence
typically the most difficult phase seems to be the inverse Z transform so it can be tricky, but we have a table, and our teacher taught us to either use

z-table
long division
partial fraction expansion and z-table

Firstly, Z-transform both sides

$Z(y[n+1])-0.5Z(y[n])=35*Z(u[n])$

according to the z-table we receive

$z*Y(z)-z*y[0]-0.5*Y(z)=35*\frac{z}{z-1}$

combine  terms and re-factor and emplace y[0]=20

$Y(z)(z-0.5)=35*\frac{z}{z-1}+z*20$

$Y(z)=35*\frac{z}{(z-1)(z-0.5)}+20*\frac{z}{z-0.5}$

We designate R(z) for the purposes of partial fraction expansion, because it is not found in the z-table

$R(z)=\frac{z}{(z-1)(z-0.5)}$

however the rightsided term is found directly in the z-table
$Z^{-1}(20*\frac{z}{z-0.5}) = 20* (0.5^k)$

start partial fraction expansions
$R(z)/z = \frac{1}{(z-1)(z-0.5)}=\frac{a}{(z-1)}+\frac{b}{(z-0.5)}$, then we find that, $a=2,b=-2$

$R(z)=2*\frac{z}{(z-1)}-2*\frac{z}{(z-0.5)}$

$Y(z)=35*(2*\frac{z}{z-1}-2*\frac{z}{z-0.5})+20*(\frac{z}{z-0.5})$

apply inverse Z transform

$y[k]=20*(0.5)^k+70*u[k]-70*0.5^k$

$y[k]=u[k]*(70-50*0.5^k)$

tabulating some results

k=-1; y[-1]=0
k=0; y[0]= 20
k=1; y1= 45
k=2; y[2]= 115/2

With such simple recursions it's often easiest (and instructive) to actually go through the recursion and find a pattern. With $y[0]=20$ and using $c=35$ we get

$$\begin{align}y[1]&=c+\frac{y[0]}{2}\\y[2]&=c+\frac{c}{2}+\frac{y[0]}{4}\\y[3]&=c+\frac{c}{2}+\frac{c}{4}+\frac{y[0]}{8}\\\vdots\end{align}$$

So we can conclude

$$\begin{align}y[n]&=c\left[1+\frac12+\frac14+\ldots+\left(\frac12\right)^{n-1}\right]+y[0]\left(\frac12\right)^n\\&=c\sum_{k=0}^{n-1}\left(\frac12\right)^k+y[0]\left(\frac12\right)^n\\&=c\cdot\frac{1-\left(\frac12\right)^{n}}{1-\frac12}+y[0]\left(\frac12\right)^n\\&=2c+\big(y[0]-2c\big)\left(\frac12\right)^n,\qquad n>0\end{align}$$

which is the same expression you obtained.


How to find PSD (Power spectral density) of spatial data I have profiled a surface, measuring the height of peaks and troughs at 0.02mm intervals. I have 1501 data points, the below listed as an example:

$$\begin{align}
x &= 0\,\mu\text m  & y &= 20\,\mu\text m\\
x &= 20\,\mu\text m & y &= 15\,\mu\text m\\
x &= 40\,\mu\text m & y &= 12\,\mu\text m\\
x &= 60\,\mu\text m & y &= 10\,\mu\text m\\
\end{align}$$

How do I generate a spatial as opposed to frequency PSD plot using Matlab? i.e. what function should I use.

Thankyou

A spectrum estimation method that works well in many cases is Welch's method. The MATLAB implementation is pretty straight forward to use, even when not knowing about the math behind it. Following Marcus' recommendations, you could try:

if norm(diff(X,2)) < sqrt(eps) % is X evenly spaced ?
    Fs = 1/mean(diff(X)); % find the step size of X
    [Pyy,F] = pwelch(Y,[],[],[],Fs,'onesided','PSD'); % compute the PSD
    plot(log10(F),log10(Pyy)), xlabel('log [\mum^{-1}]'), ylabel('log[PSD[\mum^2/\mum^{-1}]]')
    axis tight, grid on
end


Useful natural "Hilbert-like" $n$-uples and $n$-fold "analytic signals If $\mathcal{H}$ denotes the Hilbert transform, the  analytic signal of a signal $x(t)$ is $$x_a(t) = x(t) +\imath \mathcal{H}(x(t))\,.$$ The real and imaginary parts form Hilbert pairs.

Are there natural (and possibly useful) extensions, with another transform $\mathcal{T}$, which possesses involution, anti-involution or even semi-involution properties, see eg Is there a name for a function whose square is an involution?

A shape can be (for instance):
 $$x_b(t) = x(t) +\jmath \mathcal{T}(x(t)) +\jmath^2 \mathcal{T}^2(x(t))$$
of "Hilbert-like" triples? And $n$-uples?
Or, if any, related to higher order algebras (quaternions, octonions). What properties can they reveal (generalized instantaneous phase, envelopes?)

I am mainly   interested in constructions that, as least in some types of signals, in at least one application, have proven more efficient that state of the art.

The generalisation of the concept of an analytic signal is not straight forward. I'm quite certain however that looking for such a generalisation with quarternions (or even octonions) will not turn out fruitful. Those generalise complex numbers primarily algebraically, attempting to preserve as much of the field structure as possible, and not so much as a smooth geometric structure.

In my experience, the most useful generalisations of the complex numbers are usually either unitary Lie groups or Clifford algebras.

That brings me to my example, which I encountered in my work ten years ago. The mathematical background is unfortunately too vast to be laid out here, but I will attempt a sketch of the problem and explain how it relates to your question.

One of the key motivations of finding an analytic signal is the ability to extract instantaneous phase and amplitude. I will take this as the guiding theme for generalisation. In the following, I will refer to the Hilbert transform as the map from a real valued signal to its analytic counterpart and not the quadrature.

Consider the signal space $L^2(A,B)$ of square integrable functions from $A$ to $B$. The usual Hilbert transform is then a continuous (with respect to the natural topology) $\mathbb{R}$-linear shift invariant map $$H:L^2(\mathbb{R},\mathbb{R})\to L^2(\mathbb{R},\mathbb{C})$$

so that $\cos(\omega t)\mapsto \exp(i|\omega| t)$ for $\omega \in \mathbb{R}$.

We can slightly reformulate this by replacing $\mathbb{C}$ with $\mathbb{R}^+_0 \times \mathrm{U}(1)$, where $\mathrm{U}(1)$ is the unitary group of dimension 1. To be able to square-integrate that, we define the square operator $\square^2:\mathbb{R}^+_0 \times \mathrm{U}(1)\to\mathbb{R}$ so that $(a,\Phi)^2\mapsto a^2$ and $L^2(A,\mathbb{R}^+_0 \times \mathrm{U}(1))$ is well defined.

This is of course just the polar representation of the complex numbers and the Hilbert transform is still the same transform. However, with this reformulation, there is a nearly obvious generalisation. Replacing $\mathrm{U}(1)$ by a more general unitary group allows us to keep all the structure in place and have a sensible generalisation.

If we take $\mathrm{SU}(n)$ for the general case we also require an $n$-dimensional vector space for a faithful irreducible representation, so that we will allow signals in $L^2(\mathbb{R},\mathbb{R}^n)$ and the Hilbert transform is a continuous linear shift invariant map
$$H:L^2(\mathbb{R},\mathbb{R}^n) \to L^2(\mathbb{R},\mathbb{R}^+_0\times \mathrm{SU}(n) )$$
The actual map however is not obvious. It is easier to first generalise how it will be undone. The original Hilbert transform can be undone by taking the real part of the analytic result. For the generalisation, we can simply apply the group element to a unit vector, multiply the result with the amplitude and take the real part. So we map $$(a \in \mathbb{R}^+_0,\Phi \in \mathrm{SU}(n)) \mapsto a\cdot\Re\left(\Phi(\mathbf{e}_1)\right)$$ with a canonical basis vector $\mathbf{e}_1\in\mathbb{C}^n$. It is easily verified that this reduces to the standard case for $\mathrm{U}(1)$ and $\mathbf{e}_1=(1)\in\mathbb{C}^1$.

The preimages of the orbits of $\mathbf{e}_1$ generated by any generator $g\in\mathfrak{su}(n)$ of $\mathrm{SU}(n)$ under $H$ are therefore $$\Re\left( \exp(igt)\mathbf{e}_1 \right)=\frac{1}{2}\left( \exp(igt)+\exp(-igt)  \right)\mathbf{e}_1$$
where we used $g^\dagger = g$.
This suggests that $H$ should map
$$H:\frac{1}{2}\left( \exp(igt)+\exp(-igt)  \right)\mathbf{e}_1 \mapsto \left(1,\exp(igt)\mathbf{e}_1\right)$$
And here appears a problem. The map is supposed to work for all generators $g$, so also specifically for $-g$. If we plug in $-g$, the left hand side remains unchanged, but the right hand side is mapped to its complex conjugate. Hence, the map is not well defined. It is this very problem that makes using the magnitude of $\omega$ mandatory in the original Hilbert transform map. The generators in the Lie algebra form a vector space, and there is no canonical way to map all generators to a half space.

We have to equip the Lie algebra with an additional structural element, and that is a continuous map $$|\cdot|:\mathfrak{su}(n)\to\mathfrak{su}(n)$$ with the property $|g|\in \{g,-g \}$ and $|g|=|-g|$ for all $g$.

With this we have the generalisation of the Hilbert transform to unitary groups, given by the map
$$H:\frac{1}{2}\left( \exp(igt)+\exp(-igt)  \right)\mathbf{e}_1 \mapsto \left(1,\exp(i|g|t)\mathbf{e}_1\right)$$

The geometric interpretation of this transform is that a vector valued signal is decomposed into a real non-negative amplitude signal and a rotation in a complex vector space. The rotation is so the corresponding point on the Lie manifold traces a path whose time-oriented tangent is always in the same half-space, generalising the concept of positive frequency.

This does not exactly answer your original question though. You can however embed a lower dimensional signal into $\mathbb{R}^n$ and then use this generalisation to find an $n$-dimensional Hilbert transform with any number of generalised components that you desire.


Getting the power of a signal from its Fourier transform? I have a non-periodic signal that contains the sinc function in the time domain and so it is a bit difficult to calculate its power (because of the integral) through:
$$ {P_x} = \lim_{T \to \infty} \frac{1}{T} \int_{-T/2}^{T/2} |x(t)|^2 dt $$
This is my signal: $x(t)=A\cos(\omega_1t)+B\,\textrm{sinc}(\omega_2t)$. I wonder whether one can calculate the power in the frequency domain, since (Rayleigh's Theorem):
$$ \int_{-\infty}^{\infty} |x(t)|^2 dt = \int_{-\infty}^{\infty} |X(f)|^2 df $$

Note that due to the periodic component in the signal ($A\cos(\omega_1t))$, the integral
$$\int_{-\infty}^{\infty}|x(t)|^2dt=\int_{-\infty}^{\infty}|X(f)|^2df\tag{1}$$
doesn't exist, i.e., the energy of the given signal is infinite. The signal is a so-called power signal, having finite power and infinite energy.
What you can do is compute its power using the equation given in your question:
$$P_x=\lim_{T\to\infty}\frac{1}{T}\int_{-T/2}^{T/2}|x(t)|^2dt\tag{2}$$
but you can't use Parseval's theorem to compute that limit.
However, you can use Parseval's theorem indirectly to show that the contributions to the limit $(2)$ from all components of $|x(t)|^2$ containing the $\textrm{sinc}$ function go to zero as $T\to\infty$, because all those components have finite energy, hence zero power. So the power of the given signal just equals the power of $A\cos(\omega_1t)$, and I'm sure you know how to compute that power.


minimal signal window size (length) In order to detect a frequency component of 1 Hz with a DFT - do I need a signal window (signal in time domain) of at least 1 second length (in order to have at least 1 cycle of the desired minimal frequency I want to detect)?
Or is it possible to evaluate the spectrum of low frequencies with short signal windows (in seconds - independent of sampling frequency)?

If you want to detect a frequency with a DFT, yes, you'll need an integer number of signal periods in your observation. The smallest useful integer is 1.
On the other hand, why use the DFT to detect a single frequency? Aside from the obvious "use Goertzel if you only care about a single frequency", which would only lead you to a filter of the same length as said DFT, you could use one of many parametric and non-parametric frequency estimators.
A simple low- or bandpass-filter might totally do it, depending on your application!
Other estimators for when you know where you want to evaluate the spectrum include the MUSIC algorithm (which can be really low complexity if you know there's exactly one tone in your signal, and you know your SNR is OK), or very simple fast segmented autocorrelation methods.


Why Does FISTA Algorithm Not Work for Signed Signals? Using the FISTA Algorithm for compressive sensing from https://github.com/tiepvupsu/FISTA, I created the matlab example below.

I created 2 sparse signals x_signed and x_pos, where the latter only contains positive values.
Using opts.pos = false as input to the FISTA, I assume, that the algorithm also works for signed input signals.
As visible, the FISTA algorithm is able to perfectly reconstruct x_pos, but not x_signed.

s = rng(5);
close all;
nrbins = 5;
N = 120;
m = 24;
%create input signals
x_pos = [round(10*rand(nrbins,1));zeros(N-nrbins,1)];
x_signed = [10*randn(nrbins,1);zeros(N-nrbins,1)];

%random shuffle
x_signed = x_signed(randperm(size(x_signed,1)));
x_pos = x_pos(randperm(size(x_pos,1)));

%Compressive Measurements on positive signal
A = 0.5*(sign(randn(m,N))+ones(m,N));
y = A*x_pos;

opts.pos = true;
opts.lambda = 0.01;
opts.tol = 1e-14;
opts.max_iter = 1000;
xrec_pos = fista_lasso(y,A,[],opts);

%Compressive Measurements on signed signal

A = 0.5*(sign(randn(m,N))+ones(m,N));
ysigned = A*x_signed;

opts.pos = false;
opts.lambda = 0.01;
opts.tol = 1e-14;
opts.max_iter = 1000;
xrec_signed = fista_lasso(ysigned,A,[],opts);

%plot positive
figure;
scatter(1:size(x_pos,1),x_pos);
hold on;
plot(xrec_pos);

legend('x','xrec\_pos');

%plot signed
figure;
scatter(1:size(x_signed,1),x_signed);
hold on;
plot(xrec_signed);

legend('x','xrec\_signed');

How can I reconstruct a signed signal using the FISTA algorithm?

For what reason does the option "pos" exist? Does it have performance implications?

I will try explaining why would someone use such an option as defining the solution as positive.
The Proximal Gradient Method which is the hurt of the FISTA algorithm is basically a generalization of Projected (Sub) Gradient Descent Method.
As such, if one knows that the solution lies in a Convex Set and one knows the projection onto that set one could utilize it inside the method by a projection step.

Hence, if you know your solution is indeed in the convex set of the positive orthant you should use that information by projection the solution onto this set.
The projection is fairly simple:

$$ {P}_{ \mathbb{R}_{+}^{n} } \left( y \right) = \arg \min_{x \in \mathbb{R}_{+}^{n}} {\left\| x - y \right\|}_{2}^{2} = \max \left( y, 0 \right) $$

This is actually implemented in the package you linked to in the Proj files.
For instant, in the latest version of proj_l1.m you can see the code:

if opts.pos
    X = max(0, U - lambda);
else
    X = max(0, U - lambda) + min(0, U + lambda);
end

Which implements the projection onto the $ {L}_{1} $ norm ball and the positive orthant as written above.


How does the Magnitude of an FFT relate back to the Time Series wave it was created from? Alright, to start out I am a mechanical engineer with very little circuits and digital signal processing background. So please excuse my ignorance as this is far out of my usual territory.

When looking at the magnitude of an FFT of a signal, I understand that the magnitude is in the units of the time signal (in my case it is Volts). I also understand that you need to divide the FFT results by the transform length in order to normalize.

What I don't completely understand is how the magnitude of the FFT result relates back to the time signal. From what (I think) I understand it should be related to the amplitude of the time signal. From my MATLAB code that I made (posed below) it seems that it returns 1/4 of the amplitude of the signal. Is this correct? What is the physical meaning of it?

And then on a related note, if it is related to the amplitude of a time signal, what is the sampling rate needed to get a decent measure of the amplitude. I know the highest resolvable frequency for an FFT is the nyquist frequency, but I haven't been able to find anything about accurate measurements of a signal's amplitude.

Thanks in advance for any and all help on this.

MATLAB Code:

f_s = 20;              % sampling (hz)

f_w = 2;                % frequency of wave (hz)
A = 5;                  % amplitude

d_t = 1/f_s;
t = 0:d_t:2-d_t;
n = length(t);

% Creating the time signal
time_signal = A*sin(f_w*2*pi*t);

% Generating and applying a hanning window on the time signal
window = hann(n)';
windowed_time_signal = time_signal.*window;

% Taking FFT of windowed time signal
raw_fft = fft(windowed_time_signal);

% Frequency resolution of FFT
d_f = f_s/n;

% Frequency spectrum for full FFT
f_spectrum_full = (0:n-1)*d_f;

% Only first half of FFT contains useful information
half_fft = raw_fft(1:(n/2));

% Frequency scale for half FFT
f_spectrum_half = (0:n/2-1)*d_f;

% Magnitude of the half FFT
magnitude = abs(half_fft);

% Normalized version of the magnitude
normalized = magnitude/(n);

figure (1)

subplot(2,1,1)
plot(t,time_signal)
title('time wave')

subplot(2,1,2)
plot(t,windowed_time_signal)
title('windowed time wave')

figure (2)

subplot(3,1,1)
plot(f_spectrum_half,magnitude)
title('magnitude')

subplot(3,1,2)
plot(f_spectrum_half,normalized)
title('magnitude')

You can calculate the coherent or DC gain of a window with:

coherent_gain = sum(window)/n;

with a Hann window, and n = 40, this comes to 0.4875, as hotpaw2 mentioned.

You can therefore eliminate the effect of the coherent gain of the window on your spectral estimate by replacing windowed_time_signal = time_signal.*window; with windowed_time_signal = time_signal.*window/coherent_gain;.

If your main priority is very accurate amplitude - and you're unconcerned with spectral leakage as Sebastian mentioned, then the flattop window
is the way to go. Here's the matlab docs.

Just replace window = hann(n)'; with window = flattopwin(n)';.

Also, if you have the ability to collect more samples, that will likely improve the spectral estimate from your FFT (assuming the underlying process in your application is somewhat stationary/ergodic). You can also take FFTs of successive (possibly overlapping) slices and average them.

Furthermore, as hotpaw2 mentioned, selecting a sample rate depends very strongly on your application, and how many samples you intend to collect. That said, I've often used fs = 2.56 * bw, with good results when it comes to spectral analysis or digital filtering.